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I.  INTRODUCTION 


Lasers  provide  a  readily  available  source  of  coherent 
optical  radiation  for  today’s  technologies.  The  output  of 
most  lasers  exhibit  a  spatial  amplitude  distribution  which  is 
Gaussian.  It  is  possible  to  spatially  filter  such  a  beam  to 
produce  an  alternatively  shaped  beam.  Such  a  variation  may 
exhibit  a  circular  or  square  uniform  cross-section  and  either 
of  these  could  have  an  arbitrary  spatial  weighting 
distribution.  The  utility  of  such  filtering  is  unknown  unless 
the  diffracted  field  distribution  can  be  predicted  at  any 
given  distance. 

Classical  scalar  diffraction  theory  provides  a  means  for 
approximating  such  a  solution.  Scalar  theory  has  been  shown 
experimentally  to  yield  accurate  results  when  the  diffracting 
aperture  is  large  compared  to  the  wavelength  and  the  field 
locations  are  not  too  close  to  the  aperture.  Classical 
diffraction  analysis  is  a  computationally  cumbersome  method, 
however,  due  to  the  required  solution  of  multi-dimensional 
field  integrals.  For  continuous  waves,  an  alternative  method 
[Ref.  1]  expresses  the  solution  using  the  theory  of  linear 
time-invariant  systems.  If  the  multi-dimensional  Fourier 
transform  of  the  complex  field  distribution  across  any  plane 
is  taken,  the  spatial  Fourier  components  can  be  identified  as 


plane  waves  travelling  in  different  directions.  The  field 
amplitude  at  any  other  point  is  the  sum  of  each  of  these 
contributing  waves,  if  the  plane  wave's  phase  shift  during 
travel  is  accounted  for.  Thus,  the  propagation  phenomenon  may 
be  regarded  as  a  linear  space-invariant  system  characterized 
by  a  specific  transfer  function.  This  transfer  function 
behaves  as  a  linear  dispersive  filter  with  a  finite  spatial 
bandwidth.  Development  of  devices  which  generate  ultra-short 
optical  pulses  has  encouraged  the  extension  of  this  "Fourier 
optics"  concept  to  transient  waves. 

Previous  efforts  in  applying  a  transfer  function  approach 
to  transient  scalar  wave  propagation  have  focused  primarily  on 
acoustic  applications  2,3].  The  study  of  optical 
diffraction  effects  utilizing  this  mathematical  approach 
p-oves  to  be  a  natural  extension  to  that  of  acoustics. 
Indeed,  the  ultrasonic  solution  is  found  to  be  a  subset  of  the 
optical  solution. 

A  computationally  efficient  method  to  model  such 
ultrasonic  propagation  has  been  developed  by  Guyomar  and 
Powers  [Refs.  2,3].  Relying  upon  linear  systems  theory  and 
using  Fourier  transforms  to  transition  between  the  space  and 
spatial  transform  domains,  the  solution  methodology  is  an 
application  of  the  general  theory  discussed  above.  It  is,  as 
well,  quite  suitable  for  efficient  computer  implementation 
since  it  uses  Fourier  transforms.  The  mathematical  derivation 
for  the  optics  case  in  lossless,  homogeneous  material  was 
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developed  by  Guyoinar  [Ref.  4]  and  is  presented  in  detail  in 
the  followirq  chapter. 

Computer  models  of  transient  acoustic  wave  propagation 
have  previously  been  written  for  both  mainframe  and  micro¬ 
computer  applications.  An  example  of  the  latter  by  Merrill 
[Ref.  5]  was  an  early  effort  to  investigate  the  feasibility  of 
generating  solutions  on  microcomputers  within  a  "realistic" 
length  of  time  (defined  arbitrarily  and,  reasonably,  to  be 
thirty  minutes) .  An  initial  attempt  was  made  using  a 
commercially  available  program  named  MATLAB.  This  is  an 
array-oriented  program  containing  built-in  functions,  such  as 
the  Fast  Fourier  Transform  (FFT)  and  Bessel  function 
calculations.  The  capability  to  incorporate  user-developed 
functions  and  subroutines  also  exists.  (Additional  details  of 
MATLAB  will  be  presented  in  the  third  chapter  of  this  paper.) 
However,  because  MATLAB  memory  utilization  was  limited  to  64 OK 
at  that  time,  Merrill  was  forced  to  abandon  the  program  and  to 
use  Microsoft  Fortran  Version  3.31  to  develop  his  model.  The 
desired  maximum  run  time  of  thirty  minutes  was  achieved. 
[Ref.  5] 

Advances  in  the  capabilities  of  microcomputers,  including 
increased  speed  and  enlarged  RAM  capacity,  justified  a  follow- 
on  effort  to  develop  the  appropriate  code  in  MATLAB  to  be  run 
on  upgraded  processors,  such  as  the  Intel  386/486  series. 
Such  a  program  would  provide  a  convenient  means  of  researching 
optical  diffraction  field  characteristics  for  any  excitation 
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waveform.  The  specific  goal  of  this  thesis  was  thus  defined — 
to  produce  a  MATLAB  program  modelling  transient  scalar  optical 
wave  propagation.  As  no  similar  published  results  are  known 
to  the  author,  test  and  verification  of  the  output  results 
were  accomplished  by  an  informal  comparison  with  available 
unpublished  data  for  the  circular-  and  square-shaped  spatial 
excitation  waveforms  [Ref.  4].  In  addition  to  these  two 
excitation  functions,  programs  to  generate  a  spatially 
circular  waveform  with  a  selection  of  either  a  Gaussian  or 
Bessel  amplitude  distribution  were  to  be  developed  for  use  in 
further  research. 

The  following  chapter  describes  the  wave  diffraction 
problem  including  the  geometry  involved,  a  discussion  of  the 
application  of  linear  systems  theory,  and  the  mathematical 
derivation  of  the  field  solution  utilizing  the  Fourier 
approach.  Chapter  III  commences  with  an  overview  of  the 
MATLAB  language  with  particular  attention  paid  to  the 
peculiarities  in  this  application.  A  functional  explanation 
of  the  program  follows  within  the  same  chapter.  (Appendices 
C  and  D  contain  a  detailed  e:5)lanation  of  the  code.)  Chapter 
IV  illustrates  the  final  program  outputs  for  one  set  of 
defining  parameters.  Chapter  V  summarizes  the  effort  and 
discusses  areas  for  future  research.  Throughout  these 
chapters,  many  figures  are  utilized  to  enhance  the  textual 
descriptions.  These  figures  can  be  found  at  the  end  of  each 
section  within  the  appropriate  chapter  of  the  thesis. 
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Appendices  contain  further  figures,  copies  of  the  elements  of 
the  source  code,  and,  as  previously  mentioned,  a  detailed 
explanation  of  the  code's  formulation  and  operation.  Lastly, 
the  code  is  included  for  each  of  the  four  excitation 
functions. 
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II.  PROBLEM  DESCRIPTION 


The  general  solution  and  explanation  detailed  here  is 
derived  from  that  presented  by  Guyomar  and  Powers  in  Refs.  2 
and  3.  They,  in  turn,  point  out  that  the  approach  is  based  on 
the  spatial  impulse  response  introduced  by  Stepanishen  [Refs. 
6, 7, 8, 9]  to  describe  acoustic  propagation  and  reviewed  by 
Harris  [Ref.  10].  All  of  these  references  were  concerned  with 
acoustic  propagation,  but  the  optical  field  also  may  be 
expressed  as  a  temporal  convolution  of  the  time  excitation 
with  the  spatial  impulse  response.  Guyomar  and  Powers'  view 
differs  from  Stepanishen ' s  work  in  that  linear  systems  theory 
is  used  to  point  out  the  importance  of  the  total  impulse 
response  (and  its  ecjuivalence  to  the  Green's  function) .  Also, 
the  expressions  for  the  spatial  impulse  response  functions  are 
found  in  the  spatial  transform  domain.  In  this  domain, 
propagation  of  the  wave  will  be  seen  to  be  the  application  of 
two  time-varying  spatial  filters  to  the  spatial  spectrum  of 
the  input  wave. 

A.  GEOMETRY 

The  problem  utilizes  the  geometry  illustrated  in  Fig.  1. 
Given  a  separable  source  of  excitation  in  the  z=0  plane. 
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the  irradiance,  or  strength  of  the  optical  field,  at  the 
observation  plane  located  a  distance  z  from  the  source  is  to 
be  determined.  Note  that  the  excitation  is  separable  in  space 
and  time;  it  is  not  necessarily  separable  in  x  and  y. 
Propagation  through  a  lossless,  linear,  and  homogeneous  medium 
is  assumed. 


Figure  1.  Source  and  observation  planes. 


B.  LINEAR  SYSTEMS  SOLUTION  APPROACH 

It  will  be  shown  in  the  next  section  [Ref.  4]  that  the 
field  strength  for  this  separable  source  can  be  expressed  by 


(2) 


Here,  the  *  indicates  convolution  over  the  variable  appearing 
directly  below  it  and  g(x,y,2,t)  is  the  impulse  response,  or 
Green's  function,  that  both  solves  the  wave  equation  and 
satisfies  the  appropriate  boundary  conditions.  Figure  2  shows 
the  relationships  involved. 


Factoring  out  T(t)  and  its  associated  temporal 
convolution,  <^(x,y,z,t)  is  also  given  [Refs.  2,3]  as 


4)(x,y,2,t)  =  T(t)*hix,y,z,t) 

where  h(x,y,2,t)  is  the  "spatial  impulse  response"  and  is 
equal  to 


h{x,y,z,  t)  =  -s{x,y) 


dg(x,y,z,  t) 


where  g(x,y,z,t)  is  the  Green's  function  (or  impulse 
response) .  In  a  linear  system,  the  relationship  between  the 
spatial  impulse  response  h(x,y,z,t)  and  the  Green's  function 
g(x,y,z,t)  is  illustrated  by  the  block  diagram  of  Fig.  3. 
Using  the  two-dimensional  spatial  Fourier  transform  of 
Equation  4 ,  and 


Rif^.fy.z.t)  =g(4,fpFj|||  , 
0(x,y,z,t)  can  be  written  as 


(5) 
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(6) 


<|)  {x,y,  z,  t)  =  T(  t)  I 

where  the  tilde  notation  indicates  the  transform  of  the 
spatial  function. 

In  the  spatial  frequency  domain,  the  transform  of  the 
spatial  impulse  response  is 

The  two-dimensional  spatial  convolutions  found  in  Equations  2 
and  4  do  not  lend  themselves  readily  to  computer  simulation, 
due  to  the  integrations  involved.  However,  as  suggested  by 
Equations  6  and  7,  this  difficulty  can  be  overcome  by  taking 
the  appropriate  two-dimensional  spatial  transforms  using  Fast 
Fourier  Transform  (FFT)  techniques. 

Viewing  Equation  7  once  more,  s  is  the  angular  spectrum  of 
the  amplitude  distribution  of  the  input  function;  dq/dz  may  be 
thought  of  as  the  "propagation  transfer  function"  related  to 
the  medium  of  propagation.  We  will  find  that  this  propagation 
transfer  function  behaves  as  a  time-varying  spatial  filter 
that  increasingly  attenuates  the  higher  spatial  frequencies  of 
the  source  as  time  increases  [Refs.  2,3]. 

Figure  4  illustrates  the  concepts  discussed  thus  far  in 
the  time  domain.  To  summarize,  the  general  solution  technique 
begins  with  the  wave  equation  for  lossless  media  and  requires 
determining  the  Green's  function  (or  the  two  dimensional 
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spatial  transform  of  the  Green's  function)  which  solves  the 
wave  equation  and  satisfies  the  boundary  conditions. 

The  inverse  transform  of  the  product  of  the  spatial 
transform  of  the  derivative  of  the  Green's  function  and  the 
spatial  transform  of  the  input  function  will  yield  the 
solution  to  the  optical  field  values  at  the  observation  plane. 
Equation  6  alone,  or  the  combination  of  Equations  3  and  7, 
shows  this  notion  in  mathematical  form. 
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Figure  3.  Block  diagram  of  spatial  impulse  response. 


s(x,y)T(t)  Propagation 

- ►  & 

boundary  conditions 


<|)(x.y,z)  =T(t)*h(x.y,z.t) 

=-s(x,y)T(,)r,.*^^ 


Figure  4.  Block  diagram  of  general  solution. 


C.  MATHEMATICAL  DERIVATION 

The  primary  source  for  the  following  discussion  is 
Guyomar's  work  in  Ref.  4.  Some  few  additions  to  the 
mathematics  have  been  added  and  the  explanatory  text  expanded. 

The  appropriate  wave  equation,  again  assuming  a  lossless 
propagation  medium,  is  the  Helmholtz  equation  given  by 


V24) 


1  yd) 

c2  dt^ 


=  0. 


(8) 


The  fields  described  by  this  wave  equation  are  given  by  the 
radiation  integral.  However,  assuming  a  planar  aperture,  the 
field  0(x,y,z,t)  may  be  written  [Ref.  4] 


^{x,y,z,  t) 


di^o{x,y,Q,  t) 
dn 

-  4>o(-^/y/  0,  t) 


^  ^  ^ 

^9{x,y,z,t) 
X  y  t 

*  *  *  Bffix,y,z,  t) 
X  y  t  dn 


(9) 


where  0p(x,y,O,t)  is  the  excitation  function.  Its  scalar 
amplitude  distribution  is  known  on  the  source  plane.  The 
Green's  function  satisfying  the  wave  equation  and  the  boundary 
conditions  is  g(x,y,z,t)  and  the  partial  derivative  with 
respect  to  n  represents  the  normal  derivative  of  the  desired 
function  for  planar  geometry.  For  an  observation  plane 
perpendr.cu"'  to  the  z-axis,  the  normal  derivative  will  become 
tb^  derivative  A’itb.  respect  to  z. 

in  optical  applications,  the  value  of  the  field  on  the 
surface  0^(x,y,O,t)  is  known.  Therefore,  it  becomes  desirable 


12 


to  eliminate  the  first  term  of  Equation  8  and  to  work  solely 
with  the  second  term.  This  can  be  done  by  using  the  Green's 
function  which  satisfies  the  condition  that 


s\z.o  =  0 . 


(10) 


The  Green's  function  that  meets  this  condition  is 


_  +  {z-  z,y 

2 

^  ^]r^  +  (z  +  Zo)2 

_ ^ _ 

LI _ 2 _ 

+  (Z  -  Zq)^  y/l^  +  {Z  +  Zq)^ 

c  c 


(11) 


The  subscript  associated  with  the  g  Green's  function 
symbol  simply  indicates  the  sign  between  the  two  terms,  where, 
if  given  different  boundary  conditions,  a  "+"  would  be  more 
appropriate.  Note,  also,  the  conventional  use  of  the  letter 
"c"  to  indicate  the  speed  of  propagation.  On  the  surface  of 
the  source  plane  at  z=0,  this  Green’s  function  has  the 
following  properties: 


agr_  _  3gr_  (13) 

dn  ■  dz 

2z  b{t  -  R/c)  _  2z  b'{t  -  R/c)  (14) 

CR^ 

where 

Utilizing  the  Green's  function  g_,  the  field  can  be 
written  as 
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R  =  sfl^  + 


Ix^  +  +  z^ 


A/  4.  /  i-v  *  *  *  dg.{x,y,z,t) 

4^  (X/ z'/ 1)  =  ~  ^Q{,x,y f  Z/ 1)  _ - 

"  X  y  t  dz 

Substituting  Equation  13  for  dg_/dz  results  in 


4»(x,y,z,  t)  =  (j)o(x^y,0,  t)  *  *  *  2z  6  (t  -  R/c) 

X  y  t  r3 

^6  (X  V  0  t)*  *  t'{t- R/c) 


Si. ice  (f*g')  =  (f*g)  '  =  (f‘*g),  the  order  of  the  derivative  in 
the  second  term  of  Equation  18  may  be  interchanged,  yielding 

^{x,y,z,t)  =  <|)oU,y,0,  t) ~  ?-J.P1 

xyt  r3 

-  fb'  {'f  r  C  t)*  *  *  2z  6  (t  -  R/c) 

4>o(x,y,  ^  — 

Previously,  h(x,y,z,t)  had  been  defined  as  the  spatial 
impulse  response  (i.e.,  the  linear  response  to  an  input  which 
is  multi-dimensional  in  space,  yet  an  impulse  in  time) . 
Recalling  an  excitation  sourct  separable  in  space  and  time  has 
been  assumed,  this  input  is  represented  as 


4>o(x,y,0,t)  =  six, y)  bit)  . 


since  the  spatial  impulse  response  is  the  convolution  of  the 
input  function  and  the  Green's  function  (refer  to  Fig.  3) ,  the 
response  may  be  expressed  as 
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*  ♦ 


(21) 


h{x,y,z,  t)  =  s{x,y) 


X  y 


+  b'{t) 


s{x,y) 


2z  b(t  -  R/c) 

*  ♦  2z  bit  -  R/c) 
X  y  cR^ 


Rather  than  attempting  to  perform  the  spatial  convolutions 
over  X  and  y,  the  product  shall  be  taken  following  the  two- 
dimensional  spatial  Fourier  transform  of  Equation  21.  This 
transform  results  in 


Rif^,  fy,  z,  t) 


_  sif^,fy)  2zJ'o(p\/c^t^  -  z^) 


+  bU  t) 


,jZ  f-2 

sjf^.fy)  2zJ'o(pVc^t^  -  z^) 


(22) 


C^t 


where  R  is  the  two-dimensional  spatial  transform  of  h,  f^^  and 
fy  are  the  spatial  frequencies,  p  is  the  radial  spatial 
frequency  (=ff7~+~f^  ),  and  the  transform  pair  [Ref.  4] 


bit  -  R/c) 
R” 


J^ipsld^t^  -  z^) 
ict) 


(23) 


has  been  utilized. 

Recognizing  that  time  convolution  with  the  time  derivative 
of  the  delta  impulse  is  identical  with  taking  the  derivative 
of  the  function  in  the  time  domain,  i.e.. 


fit)*  b'it)  =  fHt)  , 
t 

the  spatial  impulse  response  can  now  be  written  as 


(24) 
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h{x,y,z,t)  - 

^F-^{  )H{t  -  z/c)] 

2z  J  §  ( fjf/  fy)  Jo  (  Pv^C^  ~  Z^  )H{t  -  z/c) 

c2  I  t 

(25) 

The  expression  H(t-z/c)  is  the  Heaviside  step  function  which 
indicates  that  no  output  is  possible  until  the  propagation 
wave  has  travelled  the  distance  from  the  source  plane  to  the 
observation  plane.  This  function  assures  causality  at  the 
observation  plane.  Focusing  on  the  second  term  of  Equation 
25,  the  partial  derivative  w?th  respect  to  time,  being 
independent  of  the  transform  variable,  is  moved  through  the 
transform  so  that  the  spatial  impulse  response,  in  its  final 
form,  is  expressed  as 

h(x,y,z,t)  = 

s ( f^,  fy)  -  z^  )H(t  -  z/c)] 

.  2z  a  Jo(pv/c^t^‘^)ir(t-z/c) 

- i - ^ 

(26) 

Viewing  Equation  26,  the  calculation  of  the  temporal 
impulse  response  proceeds  as  follows: 

1.  The  excitation  function  <#>j,(x,y, 0,t)  is  separated  into 
its  components  s{x,y)  and  T(t),  i.e.,  <Pg{^,Y,0,t)  = 

s(x,y)T(t) . 
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2.  The  distance  z  to  the  observation  plane  is  chosen  and 
a  value  of  time  t  is  also  chosen,  where  t  is  greater  than  or 
equal  to  z/c. 

3.  Calculate  the  two-dimensional  spatial  transform  of 
s(x,y) . 

4.  To  calculate  the  first  term  on  the  right  side, 
multiply  the  transform  by  jQ(pNc^t^  -  z^) . 

5.  Perform  a  two-dimensional  inverse-transform  on  the 
product  and  multiply  the  result  by  the  constant  2z/(c^t^). 

6.  The  second  term  requires  several  operations  prior  to 
the  inverse  transform  being  performed.  First,  divide  the 
Bessel  function  by  the  time  t.  Next,  the  time  derivative  is 
taken  and  the  result  is  multiplied  by  the  two-dimensional 
spatial  transform  of  s(x,y).  At  this  point  the  inverse 
transform  is  performed  and  its  result  is  multiplied  by  the 
constant  2z/c^. 

7.  Find  the  sum  of  the  two  terms,  and  thus  the  process  is 
complete  for  a  single  moment  in  time.  Time  may  be  incremented 
as  desired,  and  the  process  repeated  from  step  3. 

Within  the  transform  domain,  the  propagation  equation  may 
be  viewed  as  being  composed  of  two  time-varying  spatial 
"filters”.  The  first  term  of  Equation  26  shows  the 
transformed  excitation  function  being  multiplied  by  a  Bessel 
function  whose  argument  varies  with  time.  Figure  5 
illustrates  two  examples  of  the  filter  at  different  times. 
The  filter  is,  in  effect,  collapsing  in  on  itself  as  the 
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argument  grows  larger.  Additional  examples  of  this  filter  are 
shown  in  Appendix  A. 

The  partial  derivative  of  the  second  term  acts  in  a 
similar  vein.  Pictorially,  as  seen  in  Fig.  6,  it  begins  to 
appear  somewhat  similar  to  that  of  the  Bessel  function  with 
the  obvious  difference  that  the  center  of  the  function  at  p=0 
is  no  longer  the  peak  value.  Additional  examples  of  this  are 
shown  in  Appendix  B. 


Figure  5.  Time-varying  Bessel  filters  at  two  values  of 
time. 
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Figure  6.  Time-varying  derivative  filters  at  two  values 
of  time. 
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III.  PR06R2^  DESCRIPTION 


This  chapter  addresses  the  formulation  of  the  I4ATIAB 
program  written  to  model  scalar  optical  wave  propagation.  No 
in-depth  knowledge  of  MATLAB  is  assumed  and  the  discussion  of 
the  program  is  as  "functional”  as  possible.  A  detailed 
explanation  of  the  code  is  found  is  Appendices  C  and  D. 
Specific  aspects  of  MATLAB  which  drove  decisions  concerning 
the  program  is  presented  in  the  first  section.  This  is 
followed  by  a  discussion  of  the  program's  modularity,  which  is 
of  significant  importance.  A  functional  discussion  of  each  of 
the  program's  two  modules  in  separate  sections  completes  the 
presentation.  Many  figures  are  incorporated  throughout  the 
chapter  and  they  are  found  at  the  end  of  each  applicable 
section. 

A.  OVERVIEW  OP  MATLAB 

MATLAB  is  an  interactive  software  package  developed  by  The 
Mathworks,  Inc.  of  Natick,  MA.  Its  specific  orientation  is 
for  scientific  and  engineering  numeric  computations,  and 
problem  solutions  are  expressed  almost  exactly  as  they  are 
written  mathematically.  As  a  result,  many  of  the  frustrations 
and  time-consuming  development  processes  of  conventional 
programming  are  avoided.  The  name  MATLAB  stands  for  matrix 


21 


laboratory,  and  its  basic  data  element  is  a  matrix  which  does 
not  require  dimensioning.  Numerical  analysis,  matrix 
computation,  signal  processing,  and  a  graphics  capability  have 
all  been  integrated.  The  user  is  also  able  to  incorporate 
self-developed  functions,  as  has  been  done  for  the  various 
excitation  functions.  [Ref.  11] 

Of  critical  importance  to  this  project  were  the  intrinsic 
functions  of  two-dimensional  FFT's  (FFT2) ,  Bessel  function  of 
the  first  kind  of  order  zero  calculations,  and  3-D  graphics. 
The  Matlab  User's  Guide  [Ref.  10]  goes  into  some  detail  in 
explaining  the  specific  algorithm  called  by  MATLAB 's  different 
functions.  It  points  out  that  when  the  row/column  dimensions 
of  the  matrix  are  a  power  of  two,  a  high  speed  radix-two  FFT 
algorithm  is  used.  When  the  dimensions  are  not  even,  a  non- 
power-of-two  algorithm  finds  the  prime  factors  of  the 
dimensions  and  computes  the  mixed-radix  discrete  Fourier 
transform.  It  cautions  that  this  latter  process  can  become 
quite  time  consuming,  particularly  as  the  size  of  the  matrices 
become  larger.  For  this  reason,  the  decision  was  made  to  work 
with  NxN  matrices,  where  N  is  an  even  number. 

MATLAB  employs  a  backwards  three-term  recurrence  equation 
to  calculate  each  value  of  the  Bessel  function  when  the  order 
is  an  integer.  This  routine  is  time  consuming,  given  the  size 
of  the  matrices  required  and  the  number  of  values  to  be  found. 
Therefore,  the  program  execution  has  been  "modularized"  in  the 
sense  that  these  Bessel  calculations  are  accomplished 
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separately,  and  the  results,  which  are  a  function  of  radial 
distance  and  time  (see  Equation  26) ,  are  stored  to  disk  to  be 
called  when  the  second  module  is  executed.  This  has  the 
advantage  of  allowing  a  different  excitation  function  to  be 
analyzed  without  recalculating  the  Bessel  values.  This  will 
be  also  be  discussed  futher  in  an  upcoming  section. 

Additionally,  MATLAB  has  no  built-in  function  to  calculate 
the  derivative  of  a  multi-dimensional  function.  To  complete 
such  an  operation  requires  approximating  the  derivative 
through  a  difference  operation.  The  approach  taken  will  be 
explained  in  detail  in  a  following  section. 

MATIiAB's  on-line  graphics  capability  allows  viewing 
outputs  immediately,  rather  than  having  to  transport  the  data 
into  a  separate  graphics  program.  Unfortunately,  MATLAB  does 
not  numerically  scale  the  axes  when  plotting  the  data  in  3-D. 
Hence,  it  is  useful  for  comparing  relative  shapes,  but  an 
alternative  graphics  program  must  be  found  if  absolute 
amplitudes  are  to  be  displayed.  All  subsequent  figures  have 
been  generated  on  a  program  called  AXUM,  produced  by 
TriMetrix,  Inc.  of  Seattle,  Wa.  AXUM  does  have  scaled  axes, 
along  with  many  other  features,  some  of  which  will  be  briefly 
described  in  Chapter  IV. 

B.  DISCUSSION  OF  PROGRAM  MODULARITY 

As  previously  pointed  out,  the  program  has  been 
"modularized"  in  an  effort  to  separate  the  time-consuming 
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calculation  of  the  Bessel  filter  from  the  remainder  of  the 
operations. 

The  first  module  has  been  named  OPTFIL  for  OPTical  FILter. 
It  could  just  as  easily  have  been  designated  BESSFIL  for  the 
time-varying  Bessel  filter  which  it  generates,  but  work  is 
also  being  done  for  an  acoustic  case,  and  the  name  provides  a 
simple  means  of  differentiating  between  the  two. 

The  second  module  completes  the  solution  of  the  optical 
field  at  the  observation  plane  as  it  is  defined  by  Equation 
26.  This  code,  named  OPTPROP  for  OPTical  PROPagation, 
provides  for  a  selection  of  a  specific  spatially-distributed 
excitation  function,  transforms  it  via  the  FFT2  procedure, 
generates  the  derivative  type  filter,  performs  the 
multiplications  with  each  of  the  filters,  inverse  transforms 
the  results,  multiplies  each  term  by  its  respective  constants, 
and,  lastly,  adds  the  two  terms  together  for  the  output. 
Thus,  the  output  over  the  entire  observation  space  is 
computed.  However,  only  data  from  the  center  row  of  the 
resultant  matrix  (i.e.,  h(0,y,z,t))  is  used  to  construct  the 
final  output  presentation.  Each  time  increment  builds  this 
center  row  information  sequentially  as  columns  of  the  output 
data  matrix.  Such  a  display  conforms  with  the  output  plots 
most  commonly  seen  in  the  acoustics  literature. 
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C.  FUNCTIONAL  EXPLANATION  OF  OPTFIL 


OPTFIL,  in  order  to  generate  the  time-varying  Bessel 
filter,  must  initially  define  the  elements  of  the  argument  of 
the  Bessel  function.  These  elements  are  "p”,  the  spatial 
radius  of  the  filter;  "c”,  the  velocity  of  propagation  (3E8 
meters/second);  "t”,  the  incremental  value  of  time;  and  ’'z”, 
the  distance  between  the  source  plane  and  the  observation 
plane. 

Of  these,  t  requires  additional  stipulations  since  it  is 
incrementing.  The  Heaviside  step  function,  H(t-(z/c)),  is 
simulated  by  commencing  time  at  (z/c) .  The  output  plot,  it 
was  decided,  would  illustrate  the  causality  at  the  observation 
plane  by  displaying  a  small  ntunber  of  initial  time  increments 
with  zero  output.  In  the  program  this  arbitrary  number  of 
increments  is  defined  as  the  variable  "Step”.  Some  maximum 
time  must  be  specified  as  well,  and  it  is  called  *'time_max". 
The  total  number  of  time  increments  between  the  values  of 
"time  zero"  and  time_max  is  expressed  as  the  variable  "M". 
The  value  of  M  includes  the  number  of  increments  specified  in 
Step,  and  so  (M-Step)  time  "slices"  are  calculated  between  z/c 
and  time_max.  An  additional  time-related  variable  is  "t_eps". 
This  was  named  for  a  time_epsilon  added  to  each  incremental 
time  slice  for  the  purpose  of  calculating  the  derivative 
spatial  filter. 

Some  variables  are  utilized  solely  within  OPTFIL;  others 
are  needed  by  both  OPTFIL  and  OPTPROP.  Those  in  the  latter 
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category  are  stored  in  a  separate  file  called  OPTBES  which  is 
saved  to  the  hard  drive  to  be  reloaded  when  OPTPROP  begins 
executing. 

The  dimension  of  the  transform  space  and  the  calculation 
of  the  variable  "rho”  are  missing  from  the  explanation  thus 
far.  The  transform  space  is  expressed  in  terms  of  an  "NxN” 
square  matrix.  Recall  that  to  facilitate  the  FFT2 
calculation,  the  matrix  was  to  be  of  even  dimensions.  Because 
MATLAB  matrix  indices  begin  with  1  rather  than  0  (the  upper 
left  entry  being  row  1,  column  1  not  row  0,  column  0  as  an 
origin  would  require) ,  a  matrix  of  dimension  NxN  will  have  N 
points  and  N-1  segments  in  each  row  and  column.  At  first 
glance,  the  center  of  symmetry  of  the  array  would  be  at  the 
(N-H)/2  row  and  the  (N+l)/2  coltwm.  But  the  use  of  an  even 
number  of  points  means  that  this  center  of  symmetry  does  not 
have  an  array  point  associated  with  it.  Therefore,  the  center 
of  the  transform  space,  the  variable  "NO",  was  taken  to  be  at 
the  position  ((N/2)+l,  (N/2)+l).  For  a  64x64  matrix,  then, 
the  center  is  located  at  position  (33,33)  and  spatial 
frequency  radius  rho  is  measured  radially  outward  from  that 
point.  This  forced  arrangement  where  the  defined  center  is  not 
the  geometric  center  is  illustrated  with  a  circular  function 
in  Fig.  7. 

A  Fourier  Transform  (or  the  inverse  transform)  performed 
on  such  a  function  with  a  displaced  center  of  symmetry  would 
calculate  the  amplitude  correctly  but  the  phase  would  be 
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incorrect.  This  problem  with  the  phase  calculation  can  be 
overcome  by  reorienting  the  function  into  a  corner  geometry, 
as  shown  in  Fig.  8.  Here  one  can  see  the  function  centered  at 
what  would  be  the  (0,0)  position,  given  a  quadrant  I 
perspective.  The  FFT2  and  IFFT2  algorithms  within  MATLAB 
assume  a  replicated  image  beyond  the  matrix  boundaries  so  that 
mathematically  the  symmetry  of  the  function  is  retained.  The 
transform  of  the  re-oriented  function  has  the  correct 
amplitude  and  phase. 

Figure  9  depicts  a  64x64  array  base  situated  on  the  x,y 
plane  divided  into  four  quadrants.  The  notion  of  the 
quadrants  is  important  because  the  matrix  generated  by  OPTFIL 
is  symmetric  in  rho  with  respect  to  the  center  of  symmetry. 
As  a  result,  the  calculation  for  only  one  quadrant  need  be 
done  and  flipping  the  data  appropriately  fills  the  remaining 
quadrants.  This  actually  becomes  somewhat  tricky  because  of 
the  fact  that  the  quadrants  are  unequal  in  size^  Quadrant  II 
consists  of  rows  1  through  33,  and  columns  1  through  33  for  a 
33x33  "subarray”,  while  quadrant  I  includes  rows  1  through  33, 
and  columns  33  through  64  for  a  33x32  "subarray".  Quadrant 
III  is  different  still  with  dimensions  32x33,  and  quadrant  IV 
is  the  smallest  of  the  four  as  a  32x32.  Again,  the  different 
sizes  result  from  the  designation  of  element  (33,33)  as  the 
center  of  the  transform  space. 

Calculation  of  rho  is  done  in  a  straight-forward  manner 
due  to  a  MATLAB  function  called  "meshdom".  This  constructs  a 


27 


matrix  of  an  x,y  grid  incremented  as  specified  out  to  a 
specified  distance.  The  grid  elements  in  the  program  are 
named  "rhox"  and  "rhoy”.  Rhox  and  rhoy  are  the  cartesian 
sides  of  the  radial  length  rho.  Rho  was  arbitrarily  defined 
to  be  200  frequency  units  in  length  for  NO-1  poincs,  and  a 
vector  called  ”rho_m"  was  constructed  with  NO  points  so  that 
the  largest  quadrant  could  be  properly  sized.  Rho_m  is  the 
input  value  for  the  meshdom  function  and,  thereby,  an  NOxNO 
sized  matrix  is  created.  Applying  the  Pythagorean  theorem  to 
rhox  and  rhoy,  rho  at  element  (33,33)  is  zero;  at  element 
(33,64),  it  is  200;  and,  at  elements  (33,1)  and  (1,33),  it  is 
approximately  206.45.  The  maximum  value  of  approximately 
291.96  is  found  at  (1,1),  located  in  the  largest  of  tl.e 
quadrants,  as  this  location  is  the  farthest  point  from  the 
center  at  (33,33). 

To  continue  with  the  calculation  of  the  Bessel  filter,  the 
rho  value  of  each  element  in  the  x,y  grid  whose  points  are 
defined  by  (rhox,  rhoy)  is  multiplied  by  a  constant  fc^t^-z^, 
to  form  the  argument  for  the  Bessel  function.  After  the 
Bessel  operator  has  been  invoked  for  each  element,  the  results 
are  contained  within  a  temporary  matrix  called  "TEMP".  TEMP 
is  then  flipped  to  fill  in  the  quadrants  of  the  NxN  sized 
transform  space.  The  calculated  NxN  matrix  is  called  PROP. 
A  numerical  suffix  is  added  to  PROP  to  correlr  te  the  data  with 
the  time  index,  e.g.,  PROPIA,  PROP5B.  The  additions  of  *'A" 
and  "B”  serve  to  differentiate  between  those  matrices  computed 
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using  the  values  of  time  contained  within  the  time  vector  and 
those  to  which  t_eps  has  been  added  to  the  time  values.  For 
example,  PROP22A  is  computed  using  the  twenty-second  value  of 
time  in  the  time  vector.  PROP22B  is  computed  using  that  same 
time  value  plus  t_eps. 

The  variable  which  changes  from  iteration  to  iteration  is 
time.  Recall  that  the  data  commences  with  time  z/c  and 
continues  to  timejmax.  The  output  PROP  matrix  for  each  value 
of  time  is  saved  to  disk  contained  within  a  file  named 
”p(N)x(time  index) (A/B) ” .  A  numerical  suffix  is  added  to  "p” 
to  delineate  the  size  of  the  matrix  (N)  and  to  identify  the 
time  slice  from  which  it  was  formed.  For  example,  file  names 
are  of  the  form  p64xlOA  or  pl28x45B.  In  these  two  examples, 
the  matrices  are  of  dimension  64x64  and  128x128,  respectively. 
The  former  contains  the  output  matrix  PROPlOA  and  the  latter 
contains  the  output  matrix  PROP45B. 

Examples  of  the  Bessel  filters,  taken  at  time  slices  (1) , 
(3),  (8),  and  (23)  are  illustrated  in  Figs.  10-13.  The  planar 
appearance  of  the  filter  at  time  (1)  is  due  to  the  fact  that 
the  value  of  time  is  equal  to  z/c  and  the  argument  of  the 
Bessel  function  is  zero.  In  each  of  the  remaining 
illustrations,  the  maximxm  value  of  1  is  found  at  the  center 
point  of  the  transform  space  where  the  radial  distance  'ho  is 
equal  to  zero.  The  symmetry  property  is  easily  seen,  as  is 
the  progression  of  the  filter  "collapsing"  in  on  itself  as 
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Additional 


more  and  more  peaks  continue  to  be  formed, 
examples  of  the  filter  are  contained  in  Appendix  A. 
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Figure  7.  Object  with  off-center  origin  (for  viewing). 


Figure  8.  Shifted  object  (for  Fourier  transform 
calculation) . 


AMPUTUPt 


1  3: 

J  64 

QUADRANT 

II 

QUADRANT 

1 

QUADRANT 

III 

QUADRANT 

IV 

Figure  9.  Base  array  configuration. 
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Figure  11.  Bessel  filter  at  time  slice  3, 


Figure  12.  Bessel  filter  at  time  slice  8. 
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Figure  13.  Bessel  filter  at  time  slice  23. 
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D.  FUNCTIONAL  EXPLANATION  OF  OPTPROP 


This  discussion  of  the  second  program  module,  OPTPROP, 
will  be  divided  into  three  sections.  The  first  addresses  the 
input  variables  and  functions.  The  second  section  will  be  a 
functional  walk-through  of  the  code,  with  many  illustrations 
to  amplify  the  text.  These  figures  reflect  the  calculations 
at  time  (1) ,  i.e.,  only  at  time  z/c.  With  an  understanding  of 
the  variable  names  and  the  process,  the  third  section  will 
then  provide  discussion  and  illustrations  of  the  various 
computational  steps  as  time  progesses.  The  specific  time 
slices  selected  for  analysis  are  those  for  which  the  Bessel 
filter  results  were  portrayed —  time  (1),  (3),  (3),  and  (23). 
All  figures  are  incorporated  at  the  end  of  each  sub-section. 

1.  Input  variables  and  functions 

The  initial  operation  performed  by  OPTPROP  is  to  load 
from  disk  the  variables  utilized  by  both  OPTPROP  and  OPTFIL. 
Among  other  things,  these  variables  establish  the  size  of  the 
matrix,  define  the  center  point,  and  set  the  value  of  step. 
The  entire  time  vector  is  also  passed. 

OPTPROP  provides  screen  text  to  allow  a  keyboard 
selection  of  any  of  the  four  possible  excitation  functions. 
These  are  named  Circle,  Table,  Circular  Gaussian,  and  Circular 
Bessel  in  accordance  with  their  respective  spatial  amplitude 
distribution.  The  final  two  were  developed  for  future  efforts 
in  this  study,  and  no  illustrations  nor  final  outputs  will  be 
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shown  for  them.  Once  the  input  function  has  been  chosen,  the 
program  asks  for  either  the  diameter  or  the  width,  as 
appropriate.  Default  values  are  shown  in  brackets.  These 
width  values  must  be  odd  numbers  for  the  input  algorithms  to 
operate  properly,  and  the  screen  text  includes  a  reminder  of 
that. 

The  input  generated  was  viewed  via  the  built-in  MATIAB 
graphics.  An  example  of  a  Circle  with  radius  17  in  a  64x64 
matrix  is  shown  in  Fig.  14.  This  is  the  input  used  in  all  the 
following  illustrations  in  this  chapter.  Appendix  A  and 
Appendix  B. 

On  the  MATLAB  display,  the  title  "SHFT_INPUT''  win 
also  be  included.  Any  illustration,  whether  in  the  time 
domain  or  the  transform  domain  which  has  its  peak  in  the 
center  of  the  \iewing  quadrant  is  preceded  by  a  SHFT  modifier. 
This  is  to  differentiate  that  configuration  from  the  one  in 
which  preparation  for  the  FFT2  operation  has  been  done.  The 
corner  geometry  previously  mentioned  is  shown  for  this  Circle 
input  in  Fig.  15.  It,  like  every  output  which  follows  in  the 
program,  can  also  be  viewed;  its  title,  due  to  the  corner 
geometry,  is  simply  "INPUT”.  Once  it  has  been  transformed  by 
FFT2,  it  is  then  called  "F_INPUT"  for  the  frequency  domain. 
This  is  shown  in  Fig.  16.  If  shifted  back  to  the  center 
configuration  still  within  the  frequency  domain,  the  title  is 
"FSHFT_INPUT” .  See  Fig.  17.  Only  the  shifted  results  will  be 
shown  from  this  point  on. 


36 


Figure  15 .  INPUT  ( input  in  corner  geometry) . 
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Figure  16. 
geometry) . 


F_INPUT  (transfoimed  input  in  corner 


Figure  17.  FSHFT_INPUT  (shifted  transformed  input  in 
center  geometry) . 
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2.  View  at  a  single  point  in  time 

At  this  point,  the  program  begins  the  loop  through  the 
values  of  time  to  perform  the  calculations  required  by 
Equation  26.  The  equation  is  repeated  below  for  reference. 


h{x,y,z,t)  = 
2z 


F'-*{  s{f^,fy)  -  z^  )H{t  -  z/c)} 


Qf. 


(27) 


The  A  and  B  pairs  of  PROP  matrices,  containing  the  Bessel 
filter  data,  are  loaded  from  the  disk  into  the  RAM  workspace 
by  calling  the  appropriate  files,  and  are  renamed  "vnamel"  or 
'‘vname2",  respectively.  Recall  that  the  file  names  with  an  A 
addendum  to  p(N)x(time  index),  such  as  p64xlA,  contain  the 
PROP  matrices  formed  by  the  value  of  time  contained  in  the 
time  vector  at  that  index.  Those  PROP  matrices  from  files 
appended  with  a  B  utilize  the  time  value  resulting  from  adding 
t_eps  to  the  value  contained  within  the  time  vector. 

The  first  term  of  the  equation  given  above,  called 
''Fshft_outputl"  in  the  program,  is  formed  by  the  element-by¬ 
element  product  of  vnamel  and  Fshft_input.  Fshft_input  is 
necessary  rather  than  F_input  since  the  PROP  matrix  peak  is  in 
the  center  oriented  geometry.  (One  or  the  other  of  PROP  or 
F_input  required  adjustment  so  that  the  geometry  was  alike  and 
F_input  was  arbitrarily  selected.)  Fshft_outputl  for  time  (1) 
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is  shown  in  Fig.  18  (all  figures  shown  in  this  sub-section 
reflect  time  slice  (1)).  In  this  instance,  Fshft_outputl  is 
identical  to  Fshft_input  because  the  Bessel  filter  is  planar 
with  a  value  of  1.  Figure  10  shows  that  Bessel  filter. 

Calculation  of  the  second  term  was  complicated  by  the 
fact  that  MATLAB  does  not  contain  a  function  for  the 
derivative  of  a  multi-dimensional  object  such  as  the  Bessel 
filter.  For  that  reason,  a  difference  approximation  was 
adopted.  The  equation  given  below  illustrates  the  method. 
The  index  of  the  value  of  time  is  denoted  in  the  equation  by 
"m” . 

PROP  (m)  B  _  PROP  jm)  A 
t{m)+t^pg  t{m) 

In  words,  each  PROP  matrix  of  the  PROP  A/B  ^  • 
divided  by  its  respective  time.  PROPA  is  subtracted  iiom 
PROPB  and  the  result  is  divided  by  the  time  difference  between 
the  two,  which  is  t_eps.  Admittedly,  this  is  a  difference 
approximation  to  the  derivative,  and  further  work  could  be 
done  in  arranging  a  "two-sided”  derivative,  where  plus  and 
minus  t_eps  are  applied  to  a  given  reference  PROP  matrix. 
This  portion  of  the  second  term  which  contains  the  derivative 
operation  is  called  within  the  program  Der_fil  for  DERivative 
FILter.  An  example  of  this  filter  is  shown  in  Fig.  19. 
Additional  examples  are  contained  in  Appendix  B.  The  product 
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of  Der_fil  and  Fshft_input  forms  Fshf t_output2 ,  which  is  shown 
in  Fig.  20. 

In  preparation  for  the  IFFT2  operation,  both 
Fshft_outputl  and  Fshft_output2  are  converted  to  the  corner 
geometry  and  the  inverse  transform  is  then  performed.  Each  is 
multiplied  in  the  same  coding  .‘^tep  by  its  respective  constant 
to  form  "outputl"  and  ''output2''.  ”Shft_outputl"  and 
"shft_output2"  are  formed  prior  to  the  summation  step  which 
forms  "shft_output" .  Each  of  these  shifted  versions  is  shown 
in  Figs.  21-23.  Shft_output  is  not  what  will  be  actually 
observed  since  optical  detectors  are  square  law  detectors.  To 
account  for  this,  the  program  takes  the  magnitude  of 
shft_output  to  form  ''shft_outabs'' .  This  final  output  for  time 
(1)  is  shown  in  Fig.  24. 

Of  this  output  matrix,  only  the  NO  row,  h("  t) ,  is 
saved  for  inclusion  in  the  plot  which  will  depict  the  field 
over  time.  The  information  is  saved  to  one  line  of  a  matrix 
called  "output_plot" .  The  final  form  of  the  matrix  has, 
following  "Step"  columns  of  zeros,  the  results  from  the  center 
of  each  solution  matrix  for  the  successive  values  of  time  in 
successive  columns.  Thus,  the  final  output  plot  shows  the 
field  distribution  of  only  the  center  line  of  the  observation 
plane's  aperture  over  time,  although  the  time-changing  field 
over  the  entire  aperture  has  been  calculated.  Depiction  of 
the  final  output  plot  is  reserved  until  Chapter  IV.  This 
concludes  the  second  section  of  the  explanation  of  OPTPROP, 
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Figure  20.  FSHFT_OUTPUT2  (product  of  derivative  filter 
and  transformed  input  in  transform  domain) . 


Figure  21.  SHFT_0UTPUT1  (inverse  transform  of  product  of 
Bessel  filter  and  transformed  input) . 


Figure  23.  SHFT_OUTPUT  (sum  of  SHFT_0UTPUT1  and 
SHFT_0UTPUT2 )  . 
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3 .  Time  incrementing 

It  is  important  to  understand  that  the  transient  input 
function  may  be  viewed  as  a  grid  of  point  sources,  each  of 
which  is  contributing  equally  to  the  output  field.  At  any 
position  on  the  observation  plane,  the  output  field  is 
initially  formed  by  the  source  point  closest  to  that  location. 
The  field  continues  to  be  formed  as  time  piogresses  by  those 
source  points  successively  more  distant.  Thus,  the  output  is 
the  sum  of  the  effects  of  all  the  point  sources  over  a  finite 
period  of  time. 

Examples  of  the  Bessel  filter  have  been  shown  in  Figs. 
10-13.  Recall  the  description  that  the  filter  gives  the 
appearance  of  collapsing  on  itself  as  time  progresses.  The 
derivative  filter  exhibits  a  similar  tendancy  over  time,  as 
Figs.  25-28  illustrate.  Fshft_outputl,  which  is  the  element- 
by-element  product  of  the  Bessel  filter  and  the  transform  of 
the  input  function,  becomes  more  narrow  in  the  center,  and  the 
adjacent  oscillations  smooth  out  (refer  to  Figs.  29-32) .  This 
narrowing  effect  would  lead  one  to  believe  that  the  spatial 
transform  of  Fshft_outputl  should  broaden,  and  that  is  v/hat 
will  be  seen  with  shft_outputl  in  Figs.  33-36.  Fshft__output2 , 
which  is  the  result  of  multiplying  the  derivative  filter  and 
transform  of  the  input  element-by-element,  has  an  enormous 
magnitude  (see  Figs.  37-40).  Shft_oucput2  initially  is  quite 
spiked  (see  Figs.  41-44),  but  begins  to  smooth  out  over  time 
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and  remains  at  least  an  order  of  magnitude  larger  than 
shft_outputl.  For  this  reason,  the  sum  of  the  two, 

shft_output,  is  dominated  by  shft_output2  throughout.  Figures 
45-48  illustrate  this  dominance.  Shft_outabs  is  the  magnitude 
of  shft_output  (see  Figs.  49-52) . 
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amplitude 


Figure  25.  DER_FIL  at  time  slice  1  (derivative  filter) . 


Figure  26.  DER_FIL  at  time  slice  3. 
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Figure  27.  DER_FIL  at  time  slice  8 


Figure  28.  DER_FIL  at  time  slice  23. 
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Figure  29.  FSHFT_0UTPUT1  at  time  slice  1  (product  of 
Bessel  filter  and  transformed  input) . 


Figure  30.  FSHFT_0UTPUT1  at  time  slice  3. 
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input) . 


Figure  34.  SHFT_0UTPUT1  at  time  slice  3. 
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Figure  36.  SHFT_0UTPUT1  at  time  slice  23. 

54 


Figure  37.  FSHFT_0UTPUT2  at  time  slice  1  (product  of 
derivative  filter  and  transformed  input) . 


Figure  38.  FSHFT_0UTPUT2  at  time  slice  3. 


Figure  39.  FSHFT_0UTPUT2  at  time  slice  3. 


Figure  4v/.  fSHFT_0UTPUT2  at  time  slice  23. 
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Figure  46.  SHFT_OUTPUT  at  time  slice  3. 


59 


Figure  47.  SHFT_OUTPUT  at  time  yiice  8. 


60 


figure  50.  SHFT_OUTABS  at  time  slice  3. 
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Figure  51.  SHFT_OUTABS  at  time  slice  8. 


IV.  NUMERICAL  SIMULATION 


This  chapter  addressees  the  output  of  the  program  for  both 
a  circular  and  square  excitation  function,  given  a  particular 
set  of  defining  parameters.  The  parameters  used  in  the 
simulation  are  explained  in  the  first  section.  A  brief 
discussion  of  the  execution  time  follows.  The  third  section 
presents  some  information  about  the  graphics  program  AXUM 
which  was  used  to  construct  nearly  all  of  the  figures 
presentF  ■  throughout  the  thesis.  Finally,  the  optical  field 
predictions  for  each  excitation  function  are  shown. 

A.  EXPLANATION  OF  DEFINING  PARAMETERS 

The  defining  parameters  are  those  variables  found  at  the 
beginning  of  the  OPTFIL  code  which  either  OPTFIL  or  OPTPROP 
utilize  for  the  basic  problem  setup.  Examples  of  the 
parameters  include  N,  M,  Step,  time_max,  etc. ;  each  has  been 
addressed  functionally  in  previous  chapters.  This  discussion 
will  cite  the  specific  values  used  in  the  numerical  simulation 
which  produced  the  results  shown  at  the  end  of  this  chapter. 

The  parameter  N  defines  the  number  of  points  in  the  square 
matrix  and  the  value  of  N=64  was  used  for  this  case.  As  N 
defines  the  number  of  data  points,  a  128x128  matrix  would 
provide  greater  resolution  at  a  cost  of  longer  computation 
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times.  Further  research  to  investigate  whether  the  resolution 
is  warranted  would  be  useful.  The  center  of  the  matrix,  NO, 
was  33.  Step,  the  number  of  time  slices  prior  to  z/c  for 
which  a  zero  output  would  be  recorded,  was  set  at  3.  A  value 
of  64  was  selected  for  M  whl'"  ';  resulted  in  61  time  slices  (M- 
Step)  The  distance  between  the  source  plane  and  the 
observation  plane,  z.  was  set  to  80  millimeters.  This  value 
was  determined  by  scaling  the  value  of  z  used  in  prior 
acoustic  simulations  by  the  ratio  of  propagation  speeds. 
Time_max  v/as  determined  in  a  similar  fashion,  although  a  bit 
of  trial  and  error  was  also  necessary.  Its  value  in  this 
thesis  was  0.95  nanoseconds.  Rho,  also  based  on  acoustic 
work,  was  200  spatial  frequency  units,  and  c,  the  speed  of 
light,  was  set  to  3E8  meters  per  second.  T_eps  'is  the 
smallest  time  value  between  successive  Bessel  filters  for 
which  a  measurable  difference  could  be  obtained.  Anything 
smaller  meant  that  the  subtraction  of  the  two  resulted  in  a 
zero  output.  By  such  reasoning,  t_eps  was  set  to  23.5E-25 
seconds.  Table  1  summarizes  the  defining  parameters  and  the 
repective  values  used  in  this  simulation. 

B.  EXECUTION  TIME 

Although  no  formal  goal  was  set  as  far  a  desired  maximum 
execution  time  is  concerned,  it  is  naturally  of  interest.  The 
sequence  of  OPTFIL  and  OPTPROP  for  a  64x64  matrix  and  61  time 
slices  was  run  both  on  an  Intel  386/20  MHz  machine  with  a  math 


64 


Table  I.  VALUES  OF  DEFINING  PARAMETERS. 


PARAMETER 

VALUE 

DEFINITION 

N 

64 

size  of  the  square  matrix 

M 

64 

total  number  of  time  slices 

z 

80  mm 

distance  from  source  plane 
to  observation  plane 

0.95  ns 

maximum  observation  time 

rho 

200 

spatial  radius 

c 

3e8  m/s 

velocity  of  propagation 

t_eps 

23 . 5e-25s 

difference  equation  time 
value 

Step 

3 

time  increments  prior  to  z/c 

co-processor  and  on  an  Intel  486/33  MHz.  The  386/20  required 
slightly  over  forty  minutes  to  generate  all  the  Bessel  filters 
and  ten  minutes  to  complete  OPTPROP  to  the  final  output.  The 
486/33  ran  approximately  five  times  faster  needing  slightly 
less  than  eight  minutes  for  OPTFIL  and  a  bit  over  2  minutes 
for  OPTPROP. 

C.  OVERVIEW  OF  AXUM 

An  alternative  graphics  program  was  required  to  produce  3- 
D  plots  with  scaled  vertical  axes,  since  MATLAB  3-D  mesh 
displays  do  not  incorporate  that  feature. 

AXUM  is  an  advanced  technical  graphics  and  data  analysis 
program  produced  by  TriMetrix,  Inc.  It  requires  a  minimum  of 
512K  of  memory  and  generates  3-D  contour  plots  in  a  relatively 
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straightforward  fashion  [Ref.  11].  Data  is  imported  into  AXUM 
by  first  storing  the  MATLAB  data  in  ASCII  format  from  OPTPROP. 
This  was  easily  done,  and  there  are  several  examples  of  the 
coding  steps  necessary  included  in  the  source  code  in  Appendix 

F. 

Once  the  machinations  of  data  manipulation,  axes 
definition,  scaling  and  labeling,  et  al.  were  accomplished,  a 
truly  useful  feature  was  found  in  a  menu  entitled  dit 
Graph".  It  allowed  interactive  modifications  to  be  done  which 
includes  changing  the  perspective  either  horizontally  or 
vertically.  This  capability  was  used  to  orient  the 
SHFT_OUTPUT  illustrations  of  the  previous  chapter  by  rotating 
the  axes  vertically. 

D.  EX2^FLES 

The  two  source  shapes  to  be  presented  are  the  uniform 
circular  excitation  and  the  uniform  square  excitation.  Figure 
53  illustrates  the  circular  .:ase  with  a  diameter  of  17  shown 
previously  in  Chapter  3 .  Figure  54  shows  the  calculated 
spatial  impulse  response,  given  the  defining  parameters  cited 
in  Section  A  above. 

Recall  that  this  displays  the  time  progression  of 
h(0,y,z,t)  only,  for  each  matrix  generated  per  time  slice.  As 
expected,  it  is  nearly  spatially  symmetrical,  and,  over  time, 
the  strength  fades  to  zero.  "Noise"  is  present  between  the 
outboard  "tails"  and  the  center  of  the  response.  Added 
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resolution  may  smooth  that,  in  addition  to  the  other  jagged 
edges.  Note  also  the  decided  increase  in  amplitude  as  the  two 
inboard  tails  cross.  This  overall  shape  compares  favorably 
with  that  found  by  Guyomar  [Ref.  4]. 

The  uniform  square  excitation  is  shown  in  Fig.  55.  This 
particular  one  has  a  width  of  25  on  a  64x64  base.  The  output 
(see  Fig.  56)  is  smoother  than  that  of  the  circular  piston 
with  less  signal  variations  seen  throughout.  The  magnitude  of 
the  crossing  middle  tails  is  much  less  pronounced,  yet,  there 
is  an  odd  increase  in  magnitude  in  the  vicinity  of  the 
crossing.  The  fact  that  the  crossing  occurs  at  a  later  point 
in  time  can  be  explained  by  the  wider  excitation  function. 
Again,  the  field  strength  tends  to  zero  as  time  progresses. 

This  chapter  has  presented  the  optical  field  strength 
prediction  for  a  particular  set  of  defining  parameters,  using 
the  uniform  circular  and  uniform  square  excitation  functions 
of  specified  spatial  width.  The  width  of  either  excitation 
function  can  be  easily  changed  to  study  the  associated  effect 
on  the  field.  In  addition,  two  other  excitation  functions, 
with  circular  Gaussian  and  circular  Bessel  spatial 
distributions,  are  available.  The  resultant  field  strength 
predictions  using  these  latter  two  are  left  for  future 
investigation . 
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Figure  53. 
excitation. 


Input  field  for  uniform  circular  spatial 


Figure  54.  Output  field  for  uniform  circular  spatial 
excitation. 
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Figure  56.  Output  field  for  uniform  square  spatial 
excitation. 
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V.  SUMMARY 


This  thesis  investigated  the  possibility  of  solving 
complex  simulations  of  scalar  optical  wave  propagation 
utilizing  the  MATLAB  program.  Such  a  program  was  written  in 
two  sections,  OPTFIL  and  OPTPROP,  modularized  in  the  sense 
that  the  more  time-consuming  process  runs  independently  of  the 
other.  The  results  of  the  OPTFIL  are  stored  to  disk  to  be 
recalled  by  OPTPROP,  as  OPTPROP  begins  to  run.  A  change  in 
the  excitation  function,  given  a  particular  propagation 
geometry,  does  not  require  running  both  segments  of  the 
program  but  only  the  less  time-intensive  of  the  two  (OPTPROP)  . 
This  allows  for  a  more  rapid  evaluation  of  the  field 
distribution  for  numerous  trials.  Although  no  particular  goal 
for  a  maximum  run  time  was  established,  the  results  were 
satisfactory.  This  was  particularly  so  on  an  Intel  486/33  MHz 
machine  which  ran  both  segments  of  the  program  in  ten  minutes 
and  the  OPTPROP  alone  in  slightly  over  two  minutes. 

The  systems  theory  supporting  the  concept  of  a  spatial 
impulse  response  was  presented,  along  with  the  mathematical 
derivation  of  the  solution  appropriate  for  an  optical 
application.  A  functional  discussion  of  the  program,  was  given 
in  the  body  of  the  thesis  for  those  readers  whose  interest  is 
not  in  the  specifics  of  the  code.  Appendices  C  and  D  do 
contain  such  specifics.  Many  figures  have  been  incorpcrated 
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to  augment  the  discussion  of  the  program's  operation.  The 
spatial  impulse  responses  for  both  a  circular  and  a  square 
uniform  spatial  excitation  were  computed.  The  outputs  agreed 
favorably  with  previously  computed  results  [Ref.  4].  These 
responses  have  been  plotted  in  both  2-D  and  3-D  and  a  brief 
analysis  presented. 

Future  investigation  is  open  in  several  areas.  The  size 
of  the  base  matrix  should  be  increased  to  128x128  to  provide 
for  additional  resolution.  A  "two-sided"  derivative  approach 
to  the  difference  equation  may  smooth  the  results  somewhat. 
Work  thus  far  has  been  limited  to  a  square  aperture  but 
results  with  a  circular  base  could  prove  of  interest.  A  true 
understanding  of  the  propagation  will  not  be  achieved  until  a 
means  to  show  the  results  from  the  entire  field,  not  just  the 
center  portion,  is  developed.  Achieving  this,  detailed  data 
analysis  can  progress. 
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APPENDIX  A.  EXAMPLES  OF  BESSEL  FILTERS 


This  appendix  c  intains  additional  examples  of  Bessel 
filters.  The  value  of  time  at  which  each  has  been  calculated 
is  i'Ound  in  the  figure  caption. 


Figure  57.  Time  slice  1.  Figure  58.  Time  slice  3 


Figure  61.  Time  slice  23. 


Figure  62.  Time  slice  28. 


Figure  63.  Time  slice  38. 
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APPENDIX  B.  EXAMPLES  OP  DERIVATIVE  FILTERS 
This  appendix  contains  additional  examples  of  derivative 
filters.  The  value  of  time  at  which  each  was  calculated  in 
found  in  the  figure  caption. 


Figure  69.  Time  slice 


Figure  70.  Time  slice  28 


Figure  71.  Ti’-.e  slice  38. 


Figure  72.  Time  slice  58 


APPENDIX  C.  DETAILED  EXPLANATION  OF  OPTFIL 


This  explanation  addresses  in  detail  the  first  module  of 
the  program — OPTFIL.  Its  primary  function  is  to  calculate  the 
required  Bessel  filters  for  use  by  OPTPROP  in  completing  the 
diffractive  field  calculations.  This  discussion  assumes  some 
knowledge  of  MATLAB  programming  and  a  familiarity  with  the 
information  already  presented.  Refer  to  Appendix  E  for  a  copy 
of  the  code  itself. 

Initially,  OPTFIL  clears  all  variables  in  the  work  space 
and  deletes  any  previously  generated  output  files.  It  may  b-' 
desirable  to  comment  out  this  latter  command  if  only  the  size 
of  the  transfo  .  space  matrix  is  changing  and  not  any  of  the 
defining  parameters.  In  that  manner,  for  example,  64x64  sized 
PROP  matrices  may  be  retained  in  their  respective  files  while 
a  128x128  sized  solution  set  is  developed.  Given  sufficient 
space  on  the  hard  drive,  two  separate  sets  of  data  for  a  given 
combination  of  defining  parameters  are  available  to  be  used. 
Two  sets  will  require  the  defining  parameters  to  be  saved  in 
separate  files.  This  is  explained  further  below. 

Next,  the  defining  parameters  are  given  specific  values. 
This  is  not  a  user  interactive  pr  'dure  because  of  the  time 
required  to  manually  input  nine  variables  per  run  when,  in  all 
likelihood,  there  will  be  very  few  changes  once  a  set  of 
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parameters  has  been  decided  upon.  The  definition  of  each  is 
commented  to  the  right  of  its  respective  variable. 

Any  matrix  which  will  be  generated  during  the  run  of  the 
module  is  initialized  by  defining  its  size  via  the  "zeros" 
command.  The  MATLAB  User's  Manual  [Ref.  11]  points  out  that 
this  is  a  time  saving  step.  All  are  square  matrices  except 
"rho_m"  and  "time"  which  have  been  designated  vector  matrices. 

The  "time"  vector  is  next  generated  by  the  "linspace" 
command.  It  divides  the  time  period  from  z/c  to  time_max  into 
M-Step  points.  These  are  the  values  of  time  which  will  be 
used  in  the  argument  of  the  Bessel  function. 

Now  that  all  of  the  necessary  variables  have  been  defined, 
including  those  needed  by  OPTPROP,  the  program  saves  them  to 
disk  by  the  "save"  command.  The  name  of  the  file  immediately 
follows  the  comm.and  word  and  a  space  separates  the  file  name 
from  the  variable  to  be  stored  within  it.  Note  that  this 
command  allows  multiple  variables  to  be  included  within  a 
single  file.  In  this  case,  OPTBES  is  the  file  name  and  N,  M, 
NO,  Step,  time,  c,  z,  and  t_eps  are  all  stored.  If  a 
previously  generated  set  of  Bessel  functions  with  a  different 
size  matrix  is  to  be  retained  on  the  disk,  a  different  file 
name  will  be  required.  This  is  necessary  to  provide  OPTPROP 
with  the  proper  value  of  N.  OPTPROP  will  load  the  variables 
as  OPTPROP  begins  to  run. 

Calculation  of  the  spatial  radius,  rho,  was  somewhat 
tricky  because  of  the  varying  size  of  the  matrix  quadrants. 
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The  cartesian  equivalents  to  radial  values  of  rho  are 
ge..  ated  using  the  "meshdom"  command.  This  command  is 
oriented  from  the  view  of  quadrant  I.  The  shorter  axis  in  the 
case  of  an  NxN  final  matrix  has  NO-1  points  (quadrant  I  is 
NOx(NO-l)).  The  command  "linspace”  was  again  used  to  divide 
the  length  of  rho  into  NO-1  equidistant  points  to  create  the 
vector  ''rho_m''.  However,  a  length  of  N  points  exists  in  other 
quadrants  which  requires  adding  an  additional  increment  of 
"rho_m”.  The  vector  ''rho_m"  has  this  additional  increment 
added  to  it  in  the  next  line.  ''Rho_x”  and  "rhc_Y''  are  the 
cartesian  equivalents  of  the  radial  vector  ''rho_m"  and  they 
are  generated  through  the  "meshdom”  command.  A  matrix  of 
radial  distance  values  called  "row”  (to  differentiate  it  from 
"rho")  is  created  by  applying  the  Pythagorean  theorem  to 
"rhox"  and  "rhoy".  This  is  done  once  outside  the  upcoming 
loop  rather  than  repeatedly  within  the  loop. 

The  purpose  of  the  loop  is  to  sequence  the  Bessel  argument 
through  the  values  of  time  contained  in  the  time  vector.  M- 
Step  is  the  maximum  number  of  time  values  and  so  the  loop 
progresses  from  m=l  to  m=M-Step.  The  "fprintf"  command 
provides  a  screen  display  of  the  loop  count  for  an  indication 
of  the  progress  through  the  loop.  The  first  PROP  matrices  to 
be  formed  are  those  using  the  value  of  time  within  the  time 
vector.  Recall  that  these  matrices  are  delineated  through  a 
suffix  of  "A".  The  full  argument  for  the  Bessel  funcion  is 
defined  as  "arg"  and  it  is  within  "arg"  that  the  value  of  time 
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progresses.  A  temporary  matrix  called  "temp”  holds  the  result 
of  the  Bessel  operation  on  the  "arg"  matrix. 

Taking  advantage  of  the  symmetry  involved,  the  n.  xt 
sequence  of  instructions  builds  the  full  NxN  matrix  PROP  from 
the  NOxNO  matrix  "temp”.  The  first  quadrant  of  PROP  is  formed 
initially.  Since  it  has  dimensions  of  NOx(NO-l) ,  all  of  the 
row  values  of  "temp"  are  used  but  only  the  first  NO-1  columns. 
To  form  the  second  quadrant  (NOxNO)  requires  the  entire  "temp" 
matrix  "Temp"  must  be  flipped  about  the  center  to  maintain 
the  proper  relationship  of  distance,  and  the  command  "fliplr" 
accomplishes  this.  Now  that  the  entire  top  portion  of  the 
PROP  matrix  has  been  formed,  the  lower  half  is  constructed  by 
flipping  the  upper  half  about  the  center  row.  Because  the 
bottom  half  has  one  less  row  than  the  top,  only  rows  2:N0  are 
flipped. 

The  next  segment  of  the  code  correlates  the  PROP  matrix  to 
the  time  index.  The  variaole  "vname"  is  given  the  name 
"PROP (m) A".  "Vname"  is  then  set  equal  to  the  PROP  matrix 
just  formed.  This  result  is  stored  to  disk  using  the  "save" 
command  in  a  file  named  "p(N)x(m)A".  The  variables  PROP  and 
vname  are  cleared  following  this  to  save  working  space  in  RAM. 
Note  use  of  the  preparatory  command  "eval"  in  performing  these 
last  tasks. 

This  completes  the  code  which  generates  the  "A"  PROP 
matrices;  a  similar  sequence  of  commands,  with  one  exception, 
forms  the  "B"  PROP  matrices.  That  exception  is  found  within 
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the  definition  of  "arg”.  The  time  value  adds  t_eps  to  that  of 
the  value  within  the  time  vector.  Recall  that  these  "B" 
PROP’S  are  used  to  calculate  the  derivative  filter  via  a 
difference  equation  approximation. 
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APPENDIX  D.  DETAILED  EXPLANATION  OF  OPTPROP 


This  explanation  addresses  in  detail  the  second  module  of 
the  simulation — OPTPROP.  Its  primary  function  is  to  complete 
the  field  strength  calculation  once  provided  the  appropriate 
Bessel  filters  generated  by  OPTFIL.  This  discussion  assumes 
some  knowledge  of  MATLAB  programming  and  a  familiarity  with 
the  infonr.ation  already  presented.  Refer  to  Appendix  F  for  a 
copy  of  the  code  itself. 

Initially,  OPTPROP  clears  all  variables  in  the  work  space 
and  deletes  any  files  no  longer  desired.  Next,  it  loads  from 
disk  the  defining  parameters  specified  in  OPTFIL  and  contained 
in  the  file  ”OPTBES"  (if  an  additional  file  with  different 
parameter'^  has  been  created  in  OPTFIL,  ensure  the  proper  name 
called  by  OPTPROP) .  This  latter  step  is  done  very  .-.imply 
with  the  "load”  command. 

Following  this,  the  excitation  function  is  selected 
through  interactive  screen  text.  The  value  of  N  is  first 
displayed  as  a  reminder  of  the  matrix  size  and  the  choice  of 
the  four  available  input  functions  is  presented.  The  "disp” 
command  is  used  repeatedly  to  provide  the  numerous  lines  of 
screen  text.  The  desired  excitation  function  (either  circle, 
table,  circular  Gaussian  or  circular  Bessel)  is  selected  by 
keyboard  entry  of  its  respective  excitation  function  number, 
as  shown  on  the  screen  ("enter”  must  be  struck  after  the 


81 


number  is  typed) .  A  series  of  "if”  and  "elseif"  logic 
statements  follows  which  call  the  proper  function  "m"  file. 
These  function  ”m”  files  numerically  model  the  excitation 
functions  and  copies  of  their  code  are  found  in  Appendix  G. 
All  the  function  ”m"  files  require  either  the  diameter  or  the 
width  to  be  specified  as  part  of  the  call  information.  This 
information  is  again  provided  by  the  user  interactively  via 
the  keyboard.  A  default  value  is  shown  on  the  screen  in 
brackets.  (Again  "enter”  must  be  keyed  following  the 
selection.)  This  value  of  diameter  or  width  must  be  an  odd 
number  for  the  function  "m”  file  to  run  properly  and  there  is 
a  reminder  in  the  screen  text  to  that  effect.  Input  of  an 
even  value  will  result  in  an  error  message  and  the  program 
OPTPROP  will  have  to  re-initiated.  Other  information  may  be 
also  be  required  as  appropriate  to  the  particular  function. 
For  example,  if  the  circular  Gaussian  excitation  function  is 
designated,  the  screen  text  will  ask  for  the  desired  standard 
deviation,  sigma. 

The  excitation  function  is  pictorially  displayed  on  the 
screen  by  the  "mesh"  command  and  is  titled  *  SHFT_INPUT ' .  It 
is  then  changed  into  the  corner  geometry  in  preparation  for 
the  two-dimensional  Fourier  transform  operation.  The 
"fftshift”  command  performs  the  first  action  and  the  "fft2” 
command  accomplishes  the  second.  Note  that  the  result  of  any 
particular  operation  on  the  input  function  can  be  pictorially 
displayed  by  invoking  the  "mesh”  command.  As  written  in 
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Appendix  F,  the  OPTPROP  code  has  these  subsequent  "mesh" 
commands  included,  but  commented  out.  The  transformed 
excitation  function,  "F_input",  is  shifted  back  into  the 
center  of  the  transform  space  in  preparation  of  multiplication 
with  the  Bessel  filter  (the  appropriate  PROP  matrix) .  Now 
called  "Fshft_input" ,  it  is  equal  to  s(x,y)  in  Equation  26. 

The  loop  through  the  time  vector  commences  at  this  point. 
It  begins  with  m=l  and  terminates  m=M*-Step.  The  "  f print f" 
statement  provides  a  screen  display  of  the  progress  of  the 
loop.  The  "A"  file  correlated  to  the  time  vector  index  "m"  is 
"p(N)x(m)A".  A  new  variable  named  "filenamel"  is  set  equal 
to  the  character  string  "p(N)x(m)A".  file  "p(N)x(m)B"  is,  in 
essence,  renamed  "filename2*'  in  a  similar  manner.  This 
renaming  was  done  to  simplify  the  upcoming  commands. 
"Filenamel"  and  "filename2"  are  loaded  using  •»  he  "load" 
command  and  a  preparatory  "eval"  command.  Each  file  contains 
its  respective  PROP  matrix  (PROP(m)A/B) .  New  matrices,  named 
"vnamel"  and  "vname2",  respectively,  are  created  having  the 
same  values  as  PROP(m)A/B.  This  also  was  performed  to 
simplify  upcoming  commands. 

"Fshft_outputl"  is  formed  by  the  element  by  element 
product  of  "vnamel"  and  "Fshft_input" .  It  can  be  meshed  and 
the  title  will  reflect  the  loop  sequence  number  (the  same 
number  as  the  time  vector  index) .  An  example  of  the 
conversion  from  the  .mat  file  to  ASCII  format  is  shown  here. 
There  will  be  other  upcoming  steps  which  have  this  code  for 
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the  conversion  incorporated  also,  but  it  will  not  be  pointed 
out  any  further.  ”Der_fil”  is  calculated  by  the  difference 
equation  approximation  method  discussed  in  Chapter  3 . 
"Fshft__output2'"  is  the  result  of  the  element-by-element 
product  of  ’'Der_fil"  and  *'Fshft_input" .  At  this  point 
"vname"l,  '•vname2”,  ’'Der_fil'‘,  and  the  PROPA/B  pair  are 
cleared  to  free  RAM. 

To  prepare  for  the  two-dimensional  inverse  Fourier 
transform  operation,  both  ''Fshft_outputl”  and  "Fshft_output2" 
are  moved  to  the  corner  geometry,  again  by  the  "fftshift” 
command.  This  results  in  ”F_outputl''  and  ''F_output2''  being 
formed.  Both  the  inverse  transform  operation  and 
multiplication  by  the  appropriate  constants  are  done  in  the 
next  step  to  produce  “outputl"  and  ”output2''.  Each  is  shifted 
to  the  center  geometry  prior  to  the  addition  of  the  two  being 
performed.  The  sum  of  "shft_outputl”  and  ”shft_output2"  is 
called  ''shft_output".  The  ”abs”  command  is  applied  to  take 
the  magnitude  of  ''shft_output”  and  the  result  is  designated 
''shft_outabs".  Recall  that  the  magnitude  of  "sh’^t_output’'  is 
necessary  to  properly  illustrate  the  output  of  a  square  law 
optical  detector. 

At  this  point  the  NO  or  center  row  of  "shft_outabs"  is 
saved  in  the  m+Step  column  of  a  new  matrix  called  output_plot. 
MATLAB  automatically  puts  zeros  into  the  elements  of  the  first 
••Step"  number  of  columns.  (Recall  that  this  si7?ul.ates  the 
Heaviside  step  function.)  The  code  provides  the  option  to 
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similarly  view  the  magnitude  over  time  of  '•shft_outputl''  or 
”shft_output2''  to  facilitate  the  analysis  of  either. 
''Output_plotl"  and  "output_plot2"  are  the  respective  matrices 
formed.  The  end  of  the  loop  has  now  been  reached. 

When  the  loop  has  sequenced  through  all  the  values  of  time 
in  the  time  vector  and  the  matrix  "output  jplot"  is  complete, 
the  data  is  stored  to  disk  in  a  file  named  "o(d)x(M) .  The  "o" 
signifies  output,  "d"  is  the  diameter  or  the  width  of  the 
excitation  function,  and  "M"  is  the  total  number  of  time 
slices  including  Step.  An  example  of  this  filename  is 
"ol7x64".  Note  that  there  is  no  indication  of  the  specific 
type  of  excitation  function  used,  only  the  diameter  or  width 
is  specified  in  the  filename.  Th\s  file  saved  to  disk  allows 
access  later  to  the  completed  ci,tput  whether  for  view, 
plotting,  etc. 

Last  comes  a  series  of  mesh  plots  to  rovidc  viewing 
"output_plot"  from  several  aspect  angles.  The  "subplot" 
command  sets  up  two  half-screen  views  rotated  ninety  degrees 
from  each  other.  This  display  is  kept  on  the  screen  for  five 
seconds  by  the  "pause"  comm=‘  ^  before  a  full  screen  look  is 
generated.  This  compleces  OPTPROP. 
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APP)SNDIZ  E.  SOURCE  CODE  FOR  OPTFIL 


This  appendix  contains  the  source  code  for  OPTFIL.  The 
code  is  written  using  the  MATLAB  software  package. 

OPTFIL 


%%  OPTFIL. M 

%%  This  program  generates  the  Bessel  filter. 

%%  June  6,1992 
clear; 

%!del  p*x*.mat 

N  =  £4;  %  size  of  square  array 

M  =  64;  %  number  of  time  slices 

NO  =  (N/2)+l;  %  defines  center  of  the  square  array 

Step  =3;  %  number  of  leading  zero  time  slices. 

2  =  80e-3;  %  distance  to  the  observation  plane 

time_max  =  .95e-9;  %  time  at  the  final  time  slice 

rho  =  200;  %  spatial  radius  of  the  filter 

tsqrt(rhox^2 

%  +  rhoy^2) ] 

c  =  3e8;  %  velocity  of  the  light  wave 

=  23.5e-25;  %  time  difference  for  the  differential 

%%  Initialize  matrices  to  save  processing  time 

PROP  =  zeros (N) ; 

temp  =  zeros (NO); 

arg  =  zeros (NO); 

rho_m  =  zeros (NO, 1); 

row  =  zeros (NO) ; 

time  =  zeros (M-Step, 1) ; 

%%  Ge.ierate  M-Step  time  slices  between  z/c  and  time_max. 
time  =  linspace( z/c, time_max, M-Step) ; 

%%  Save  those  variables  necessary  for  OPTPROP.m  in  a  file 
%%called  OPTBES. 

save  OPTBES  N  M  NO  Step  time  c  z  t_eps; 

%%  Generate  NO-1  values  of  rho_m  from  0  to  rho. 
rho_m  =  1 inspace (0, rho, NO-1) ; 

%%  Add  additional  increment  to  rho_m  to  compensate  for  the 
%%off-center  orientation  of  the  final  NxN  matrix 
rho_m  =  [rho_m  (rho_m(N0-l)+rho_m(2) ) ] ; 


86 


%%  Create  two  NO  x  NO  arrays  of  rho  values  for  -^unction 
%%evaluation. 

[  rhox ,  rhoy  ]  =  meshdom  ( rhojm ,  rho_in) ; 

%%  Calculate  matrix  of  radial  distance  values  outside  the  loop 
row  =  sqrt(rhox.^2  +  rhoy. ^2) ; 

for  m  =  l:M  -  Step  %%%%%%%%%%START  LOOP%%%%%%%%%% 

fprintf ('%3.0f ’ ,m) ;  %show  m  value  on  screen 

%%  Generate  PROP  matrices  with  suffix  of  ”A”  corresponding 
%%  to  the  values  of  the  time  vector 

%%  Create  an  NO  x  NO  array  of  argument  values  for  the 
%%  bessel  function 

arg  =  row  *  sqrt(  c^2* (time(m) ) ^2-z^2  ); 

%  Evaluate  J__0  at  each  argument  value;  creates  an  NO  x  NO 
%  array 

temp  =  bessein(0,arg) ; 

%  Create  PROP  matrix  containing  the  Bessel  filter  data 
PROP(1:NO,NO:N)  =  temp(l:NO,l;NO-l) ; 

PROP(1:NO,1:NO)  =  fliplr(temp) ; 

PROP(NO;N,1:N)  =  flipud(PROP(2:NO,l:N) ' 

%Correlate  the  name  of  the  variable  PROP  with  the  time 
%index;ie,  PROPIA,  PR0P2A, . . . 

vname  =  [ 'PROP' , int2str(m) ,  'A' ] ;  %set  up  name 

eval ( [vname, '=  PROP  ; ' ] ) ; 

%oave  applicable  PROP  in  a  file  named  p(N)x(m)A; .e.g. , 
%PR0P5A  in  p64x5A 

eval{[ 'save  p' ,int2str(N) , 'x' ,int2str(m) , 'A' , '  ' , vname]) ; 
eval ( [ ' clear  PROP  ' , vname] ) ; 

%Generate  PROP  matrices  with  suffix  of  "B"  corresponding 
%to  the  values  of  the  time  vector  +  t_eps 

%  Create  an  NO  x  NO  array  of  argument  values  for  the 
%  be-sel  function 

arg  =  row  *  sqrt(  c^2* (time(m)+t_ops) ^2-z^2  ); 

%  Evaluate  J_0  at  each  argument  value;  creates  an  NO  x  NO 
%  array 

temp  =  besseln(0,arg) ; 

%  Create  PROP  matrix  containing  Bessel  filter  data 
PROP(1;NO,NO:N)  =  temp(i:NO,l:NO-l) ; 

PROP(l:NO,l;NO)  =  fliplr(temp) ; 

PROP(NO:N,1:N)  =  flipud(PR0P(2 :N0, 1:N) ) ; 
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end 


%Correlate  the  name  of  the  variable  PROP  with  the  time 
% index ; ie ,  PROPIB ,  PR0P2  B , . . . 

vname  =  [ 'PROP* , int2str (m) , 'B' ] ;  %set  up  name 

eval ( [vname, ‘=PROP  ; ' ] ) ; 

%Save  applicable  PROP  in  a  file  named  p(N)x(m)B;e.g. , 
%PR0P6B  in  p64x6B 

eval(['save  p%int2str(N) , 'x*  ,int2str(m) , 'B' , '  vname]); 
eval ( [ ' clear  PROP  ’ , vname] ) ; 

%%%%%%%%%%%END  LOOP%%%%%%%%%% 
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AFFEMDIZ  F.  SOURCE  CODE  FOR  OFTFROF 


This  appendix  contains  the  source  code  for  OPTPROP. 
code  is  written  using  the  MATLAB  software  package. 
OPTPROP 


%%  OPTPROP. m  performs  transient  optical  wave  propagation 
%%  simulations.  It  uses  the  NxN  arrays  ''p(N)x(m)A/B"  to 
%%  compute  the  propagation  transfer  funct^^on. 

%%  June  2,  1992 

%%  Size  of  the  variables 

%%  NxN  —  input  functions;  M-Step  —  time  slices. 

%%  NxM — output_plot 

clear; 

idel  opt*. met 
%!del  optder?.dat 


%%  Load  the  defining  parameters  specified  in  OPTFIL.m 
load  OPTBES.MAT 


%%  Generate  the  INPUT  function;  plot  it. 
N 


’) 

') 


disp('N  is  the  width  of  the  base  for  each  function') 
disp('  ') 

disp(' Please  select  the  excitation  function') 
disp(' 
disp( ' 
disp ( ’ 

disp ( '  4  -  Circular  Bessel  ' ) 

disp( ' 
disp( ' 
disp( ' 
disp( ’ 


1  -  Circle 

2  -  Table 

3  -*  Circular  Gaussian 

4  -  Circular  Bessel 


’) 

Please  strike  "Enter"  after  selection.') 
') 

Default  values  are  in  [  ] . ‘ ) 


•) 


input_f unc  =  input ( ' Please  enter  an  excitation  function 

number  [1]  ’); 
if  isempty {input_func) 
input_func  -  1 

end 


if  input_func  ==  1, 

d  =  input ( 'Please  enter  ODD  diameter,  [13],  d  =  '); 

if  isempty (d) 
d  =  13 
end 


The 
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shft_input  =  crcle(d,N); 
elseif  input_func  ==  7., 

w  =  input (' Please  enter  ODD  width,  [25],  w  =  •); 

if  iseinpty(w) 
w  =  25 
end 

shft_input  =  table (w,N); 
d  =  w; 

elseif  input_func  =  3, 

sigma  =  input {• Please  enter  the  standard  deviation,  [12], 
sigma  =  ' ) ; 

if  isempty (sigma) 
sigma  =  12 
end 

d  =  input (' Please  enter  the  ODD  diameter,  [25],  d  =  '); 

-if  isempty  (d) 
d  =  25 
end 

shft_^input  =  crcgaus ( sigma, d, N) ; 
elseif  input_func  ==  4, 

a  =  input ( 'Please  enter  the  width  scaling  factor,  [1] ,  a=  ' )  ; 
if  isempty (a) 
a  =  1 
end 

d  =  input ( 'Please  enter  the  ODD  diameter,  [25],  d  =  '); 

if  isempty (d) 
d  =  25 
end 

shft_input  =  crcbess(a,d,N) ; 
else 

error ( ' Incorrect  Excitation  Function  Selection ' ) 
end 

mesh(shft_input) ; title ( 'SHFT_INPUT' ) ; 

%following  code  used  to  enable  print  function  and  save 
%data  to  disk 
%meta  optin 
%pause(l) 

%clg 

%vname  =  [ ' shf t_input ' ] ; 

%eval ( [ *  save  opshftin  ' , vname] ) 

%eval ( [ ' save  opshftin . dat  ' , vname , '  /ascii ' ] ) 

%%  Shift  input  quadrants  and  take  the  2-D  FFT  to  produce 
%%  F_INPUT. 

input  =  (fftshift(shft_input) ) ; 

%following  code  used  to  observe  and  save  data  to  disk 
%mesh ( input) ; title ( ' INPUT  ' ) 
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%pause(l) ;clq 
%vnaine  =  [ '  input '  ] ; 

%eval ( [ ' save  opt in  ' , vname] ) 

%eval(['save  optin.dat  vname, ’  /ascii’]) 

F_input  =  fft2 (input); 

%mesh ( F_input ) ; title ( ’ F_INPUT  ' ) 

%pause(l) ;clg 

%%  Shift  F_input  in  preparation  of  multiplication  with  PROP 
Fshft_input  =  fftshift(F_input) ; 

%mesh (Fshf t_input) ; title ( ’ FSHFT_INPUT  ' ) 

%pause ( 1 ) ; meta  opt f_in ; clg 

%%  Plot  the  absolute  value  of  the  (shifted)  transform  for 
%%  viewing  only 

%mesh(abs (Fshft_input) ) ;  title( 'ABS (FSHTFT_INPUT) ’ ) 
%pause(l) ;clg 


%%  Array-multiply  the  shifted  transfer  function  PROP  and 
%%  Fshft_input. 

disp( 'Performing  array  multiplication. ...'); 
for  m  =  1:M-Step 

%%%%%%%%%%%%  Start  loop  %%%%%%%%%%%%%%%%%%%%% 
fprintf ( ' %2 . Of , ' ,m) 
pause (1) 

filenamel  ==  [ 'p' ,int2str(N) , 'x' ,int2str(m) , 'A'  ]; 
eval ( [ ' load  ' , filenamel]  ) ; 

filename2  =  [ 'p* ,int2str(N) , 'x' ,int2str(m) , 'B'  ]; 
eval (['load  ',filename2]  ); 
eval ( [ ' vnamel  =  PROP' , int2str (m) , 'A' ,';']) ; 
eval ( [ 'vname2  =  PROP' , int2str(m) , 'B' ,';']) ; 

Fshft_outputl  =  (vnamel  .*  (Fshft_input) ) ; 

%following  code  used  to  observe  at  given  time  slice  (m) 
%and  save  data  to  disk 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


if  m  ==  46 


mesh ( Fshf t_outputl) ; 

title ( [ ' FSHFT_0UTPUT1  ( ' , int2str (m) , ' ) ' ] ) 
pause (1) 
meta  optl; 

vname  =  [ ' F_outl ' , int2str (m) ] ; 
eval ( [vname, ' =Fshf t_outputl  ; ' ] ) 
eval ( [ ' save  optf 1_' , int2str (m) , '  ' , vname] ) 
eval ( [ ' save  optf 1_' , int2str (m) , ' . dat  ' , vname , ' 
/ascii' ] ) 


end 

clg 


Der_fil =  ( (vname2/(time(m)+t_eps) )  -  (vnamel/time (m) ) ) . . . 

/t_eps ; 

%following  code  used  to  observe  at  given  time  slice  (m) 
%and  save  data  to  disk 
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if  m  ==  46 

mesh(Der_fil) ; 

title ( [ • DER_FIL  ( • , int2str (m) , • ) • ] ) 
pause (1) 
meta  optdir; 

vname  =  [ ' Der_fil ’ , 5nt2str (m) ] ; 
eval  ( [vname , '  =Der  f  ii  i* '  ] ) 
eval  ( [  'save  optdir' ,  ir.t2str  (m) , '  ' ,  vname] ) 
eval ( [ ' save  optdir • • int2str (m) , ' . dat  ' , vname , ' 
/ascii' ]) 

end 

clg 

Fshft_output2  =  Der_fil  .*(Fshft_input) ; 

%following  code  used  to  observe  at  given  time  slice  (m) 
%and  save  data  to  disk 
if  m  -=  46 

mesh (Fshf t_output2 ) 

title { [ ' FSHFT_0UTPUT2  ( ' , int2str (m) , • ) ' ] ) 
pause (1) 
meta  opt2 

vname  =  [ ' F_out2 ' , int2str (m) ] ; 
eval ( [vname , ' =Fshf t_output2  ; ' ] ) 
eval(['save  optf2_' ,int2str(m) , '  vname]) 
eval( [ 'save  optf2_' ,int2str(m) , ' .dat  ' , vname, ' 
/ascii*  3 ) 

end 

clg 

%%  Clear  unnecessary  variables  to  free  RAM 
clear  vnamel;  clear  vname2;  clear  Der_fil; 
eval ( [ 'clear  PROP' ,int2str{m) , 'A* ,';']); 
eval ( [ 'clear  PROP' ,int2str(m) , 'B' ,';']) ; 

%%  Shift  Fshft_output(l,2)  to  corner  geometry  prior  to 
%%taking  the  IFFT2 

F_outputl  =  fftshift(Fshft_outputl) ; 

mesh(F_outputl)  ;  title (' F_0PTPUT1  ') 

pause (1) ; 

clg; 

F_output2  =  f ftshift(Fshft_output2) ; 

mesh(F_output2)  ;  title (* F_0UTPUT2  '); 

pause (1) ; 

clg; 

%%  Take  IFFT  of  F_output ( 1 , 2 )  and  multiply  each  by  its 
%%  respective  constants  to  produce  the  outputs 
outputl  =  (ifft2  (F_outputl) )  *(  2*z/  ( (c^2)  *  (time(m)  ^2) )  ); 
mesh(outputl) ;  title (' OUTPUTl  '); 
pause (1) ; 
clg; 
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OUtput2  =  (if ft2 {F_output2) ) * (2*2/ (C^2) ) ? 

mesh (output 2)  ;  title (' 0UTPUT2  •); 

pause (1) ; 

clg; 

%%  Shift  output (1,2)  prior  to  summation 
shft_outputl  =  fftshift(outputl) ; 

%following  code  used  to  observe  at  given  time  slice  (m) 
%and  save  data  to  disk 
if  m  ==  46 

mesh ( shf t_output 1 )  ; 

title ( [ ' SHFT_0UTPUT1  ( • , int2str (m) , ' ) * ] ) 
pause (1) ; 
meta  optoutl 

vname  =  [ ’ outl ' , int2str (m) ] ; 
eval ( [vname , ’ =shf t_outputl  ; ' ] ) 
eval ( [ 'save  optl_' , int2str (m) , '  ' , vname] ) 
eval ( [ ' save  optl_ * , int2str (m) , ' . dat  ' , vname , ' 
/ascii' ] ) 

end 

clg; 

shft_output2  =  fftshift (output2) ; 

%  foil  owing  code  used  to  observe  at  given  time  slice  (m) 
%and  save  data  to  disk 
if  m==  46 

mesh ( shf t_output2 )  ; 

title  ( [ '  SHFT__0UTPUT2  ( ' ,  int2str  (m) , ' )  » ] ) 
pause (1) 
meta  optout2 

vname  =  [ 'out2 ' , int2str (m) ] ; 
eval ( [vname, '=shft_output2  ; ' ] ) 
eval([ 'save  opt2_' ,int2str(m) , '  ' , vname] ) 
eval ( [ ' save  opt2_' , int2str (m) , ' . dat  ' , vname , ' 
/ascii' ] ) 

end 


%%  Calculate  the  shifted  output 
shft_output  =  shft_outputl  +  shf t_output2 ; 

%following  code  used  to  observe  at  given  time  slice  (m) 
%and  save  data  to  disk 
if  m  ==  46 

mesh ( shf t_output) ; 

title ( [ ' SHFT_OUTPUT  ( ' , int2str (m) , ' ) ' ] ) 
pause (1) 
meta  optout 

vname  =  [ ' out ' , int2str (m) ] ; 
eval ( [vname , ' =shf t_output  ; ' ] ) 
eval ( [ ' save  opt3_' , int2str (m) , '  ' , vname] ) 
eval(['save  opt3_' , int2str (m) , ' .dat  ', vname,' 
/ascii' ] ) 
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%  end 

%  clg 

%%  View  the  magnitude  of  the  shifted  output 
shft_outabs  =  abs(shft_output) ; 

%following  code  used  to  observe  at  given  time  slice  'm) 
%and  save  data  to  disk 

%  if  m  ==  46 

%  mesh(shtt_output) ; 

%  title ( [ ‘MAG  SHFT_0UTPU1  ( ' , int2str (m) , ' ) * ] ) 

%  pause (1) 

%  meta  optabs 

%  vname  =  [ 'outnbs* ,int2str(m) ] ; 

eval ( [vname#  *=shft _outabs  ; ' ] ) 

eval  ( [ '  save  op'cab_  ^ ,  int2  str  (m) , '  ' ,  vname  ] ) 

%  eval ( C ' save  optab_ • , int 2  str (m) , ‘ . dat  ' , vname , ’ 

/ascii' ] ) 

%  end 

%  clg 

%%  Save  the  NO  j  if  the  m.:^nitude  of  the  shifted  ouvputl 

%%  in  the  m/2  cc  n  of  ou"  JUt  plotl. 

output__plotl (1;N,  •  »tev';  -  a;  e;''hft_outputl(NO,l:N) )  *  ; 

%%  Save  the  NO  row  the  :  xcude  of  the  shifted  output2 

%%  in  the  m/2+Step  column  ot  output_plot2 . 
output__plot2  (l:N#m->  Step)  =  abs(shft_output2  (VO,  1:N) )  * ; 

%%  Save  the  NO  row  oZ  .he  magnitude  of  the  shifted  output 
%%  in  the  m+Step  column  of  output_plot. 
output_plot(l;N,m+Step)  =  ■.>hft_outabs(NO,  1:N)  ' ; 


end 


End  loop  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


00000 

2:  ^  sr  Jg  is: 

00000 


%  Save  contents  cf  "output_plot"  as  NxM  array, 

filename  =  [ 'o' , int2str (d) , 'x' , int2str(M) ] ;  %  File:  o(d)x(M) 

eval ( [ ' save  ' , filename, '  output_plot ' ]  ) ; 

%%  Plot  the  absolute  value  of  "outputjplot" . 
subplot (121) ,  mesh ( rots O(output_plot,l) ) ; 
subplot (122) ,  mesh(rot90(output_plot,2) ) ; 
pause (5) 

%%  Plot  full  screen  view 
clg; 

mesh (rotSO (output_plot, 3) ) ; 
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APPENDIX  6.  SOURCE  CODE  FOR  EXCITATION  FUNCTIONS 


This  appendix  contains  the  source  code  for  each  of  the 
four  excitation  functions.  These  functions  are  the  uniform 
circular  and  uniform  square  functions,  and  the  circular 
Gaussiaii  and  circular  Bessel  functions. 

UNIFORM  CIRCULAR  EXCITATION  FUNCTION 


function  Y  =  crcle(d,N) 

%  CRCLE.M:  Y  =  crcle(d,N) 

%  Program  for  generating  uniform  circular  excitation  functions 
%  d  is  the  DIAMETER  of  the  circle.  (ODD  integer) 

%  N  is  the  WIDTH  of  the  square  base.  (EVEN  integer) 

%  Example;  z  =  crcle(33,64) ; 

%%%%%  JG  Upton  March,  1992 

%  Check  that  d  is  an  odd  integer 
if  rem(d,2)  <  0.1; 

error ( ' The  diameter  of  the  crcle  function  must  be  an  ODD 
integer. ' ) ; 
else; 
end; 

%  Check  that  N  is  an  even  integer 
if  rem(N,2)  “=  0.0; 

error  (’The  width  of  the  square  base  must  be  an  EVEN 
integer. ’ ) ; 
else; 
end; 

NO  =  (N/2)+l; 
r  =  d/2; 

Y  =  zeros (N) ; 
temp  =  zeros (NO-1); 

for  m  =  l;r+l; 

for  n  =  l:r+l; 

if  sqrt(  (m-1) -^2  +  (n-l)^2)  <=  r; 
temp(m,n)  =  1; 
end; 

end; 

end; 


%N0  is  the  center  location 
%r  is  the  radius 
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Y(NO:N,NO:N)  =  temp; 

Y(2:N0,N0:N)  =  flipud(temp) ; 
Y(2:N0,2:N0)  =  rot90 (temp, 2) ; 
Y(N0;N,2;N0)  =  fliplr(temp) ; 
%%%%%For  testing  purposes;  mesh(Y) 


UNIFORM  SQUARE  EXCITATION  FUNCTION 


function  Y  =  table {w,N) 

%  TABLE. M;  Y  =  table (w,N) 

%Program  for  generating  a  uniform  square  excitation  function. 
%  w  is  the  WIDTH  of  the  table.  (ODD  integer) 

%  N  is  the  WIDTH  of  the  square  base.  (EVEN  integer) 

%  Example:  z  =  table(33,64) ; 

%%%%%  JG  Upton  March,  1992. 

%  Check  that  w  is  an  odd  integer, 
if  rem(w,2)  <  0.1; 

error  ( '  The  width  of  the  table  must  be  an  ODD  integer . ' )  ? 

else; 

end; 

%  Check  that  N  is  an  even  integer 
if  rem(N,2)  “=  0.0; 

error  (’The  width  of  the  square  base  must  be  an  EVEN 
integer. ’ ) ; 
else; 
end; 

NO  =  (N/2)+l;  %  NO  is  the  center 

location 

wO  =  ceil(w/2)  ;  %  wO  is  the  mid-point  of 

the  table 

Y  =  zeros (N); 

temp  =  zeros (NO-1); 

temp(l;w0, l:w0)  =  ones(wO); 

Y(N0:N,N0:N)  =  temp; 

Y(2:N0,N0:N)  =  rot90(temp); 

Y(2:N0,2:N0)  rot90(temp,2) ; 

Y(N0;N,2:N0)  =  rot90(temp,3) ; 

%%%%%For  testing  purposes:  mesh(Y) 


CIRCULAR  GAUSSIAN  EXCITATION  FUNCTION 


function  Y  =  crcgaus (sigma, d,N) 

%  CRCGAUS. M;  Y  =  crcgaus (sigma, d,N) 
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%Prograin  for  generating  circular  Gaussian  ; Itation 
functions. 

%  sigma  is  the  STANDARD  DEVIATION  of  the  gaussian  function. 
%  d  is  the  DIAlffiTER  of  circle.  (ODD  integer) 

%  N  is  the  WIDTH  of  the  square  base.  (EVEN  integer) 

%  Example:  z  =  crcgaus(12,33,64) ; 

%%%%%  JG  Upton  March,  1992. 


iftu=0;  %mu  is  the  mean  of  the 

gaussian  function 

%  Check  that  d  is  an  odd  integer, 
if  rem(d,2)  <  0.1; 

error ( • The  diameter  of  the  circle  function  must  be  an 
ODD  integer . ' ) ; 
else; 
end; 

%  Check  that  N  is  an  even  integer, 
if  rem(N,2)  “=  0.0; 

error  ( ’  The  width  of  the  square  base  must  be  an  EVEN 
integer. ' ) ; 
else; 
end; 

NO  =  (N/2)+l; 
r  =  d/2; 

Y  =  zeros (N); 
temp  =  zeros (NO-1) ; 

for  m  =  1: (d+l)/2? 

for  n  =  l:(d+l)/2; 

X  =  sqrt((m-l)^2+(n-l)^2) ; 
if  X  <=  r; 

temp(m,n)  =  (l/(sqrt(2*pi) *sigma) ) *exp(-( (x-mu) ^2)/. . . 
(2*(sigma^2) ) )  ; 

end; 

end; 

end; 

Y(N0:N,N0;N)  =  temp; 

Y(2:N0,N0:N)  =  flipud(temp) ; 

Y(2:N0,2:N0)  =  rot90 (temp, 2) ; 

Y(N0:N,2:N0)  =  fliplr(temp) ; 

%%%%%For  testing  purposes:  mesh(Y) 


%  NO  is  center  location; 
%  r  is  the  radius 
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CIRCULAR  BESSEL  EXCITATION  FUNCTION 


function  Y  =  crcgaus (sigma, d,N) 

%  CRCGAUS.M:  y  =  crcgaus (sigma, d,N) 

%Program  for  generating  circular  Gaussian  excitation 
functions. 

%  sigma  is  the  STANDARD  DEVIATION  of  the  gaussian  function. 
%  d  is  the  DIAMETER  of  circle.  (ODD  integer) 

%  N  is  the  WIDTH  of  the  square  base.  (EVEN  integer) 

%  Example:  z  =  crcgaus (12, 33, 64 ) ; 

%%%%%  JG  Upton  March,  1992. 


mu=0;  %mu  is  the  mean  of  the 

gaussian  function 

%  Check  that  d  is  an  odd  integer, 
if  rem(d  2)  <  0.1; 

error ('The  diameter  of  the  circle  function  must  be  an 
ODD  integer.'); 
else; 
end; 

%  Check  that  N  is  an  even  integer, 
if  r6m(N,2)  *=  0.0; 

error  ('The  width  of  the  square  base  must  be  an  EVEN 
integer. ' ) ; 
else; 
end; 

NO  =  (N/2)+l; 
r  =  d/2; 

Y  =  zeros (N); 
temp  =  zeros (NO-1); 

for  m  =  1: (d+l)/2; 

for  n  =  l;(d+l)/2; 

X  =  sqrt((m-l)''2+(n-l)^2)  ; 
if  X  <=  r; 

temp(m,n)  =  (l/(sqrt(2*pi)*sigma) )*exp(-((x-mu)^2)/. . . 
(2* (sigma^2) ) ) ; 

end; 

end; 

end; 

Y(N0:N,N0;N)  =  temp; 
y(2;N0,N0:N)  =  flipud(temp) ; 
y(2;N0,2:N0)  =  rot90(temp,2) ; 
y(N0;N,2;N0)  =  fliplr(temp) ; 

%%%%%For  testing  purposes;  mesh(Y) 


%  NO  is  center  location; 
%  r  is  the  radius 
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CIRCULAR  BESSEL  EXCITATION  FUNCTION 


function  Y  =  crcbess(a,d,N) 

%  CRCBESS.M:  y  =  crcbess{a,d,N) 

%  Program  for  generating  circular  Bessel  excitation  functions. 
%  a  is  the  WIDTH  SCALING  FACTOR. 

%  d  is  the  DIAMETER  of  the  circle.  (ODD  integer) 

%  N  is  the  WIDTH  of  the  square  base.  (EVEN  integer) 

%  Example;  z  =  crcbess(l,33,64) ; 

%%%%  J.G.  Upton  March,  1992 

%  Check  that  d  is  an  odd  integer, 

if  rem(d,2)  <  0.1; 

error  ( ’  The  diameter  of  the  circle  must  be  an  ODD 
integer ' ) ; 

else; 

end; 


%  Check  that  N  is  an  even  integer, 

if  rem(N,2)  “=  0.0; 

error  ( '  The  width  of  the  square  base  must  be  an  EVEN 
integer ' ) ; 

else; 

end; 

NO  =  (N/2)+l;  %N0  is  the  center 

location 

r  =  d/2;  %r  is  the  radius  of  the 

circle 

y  =  zeros (N) ; 
temp  =  zeros (NO-1); 


for  m  =  l:r+l; 

for  n  -  l;r+l; 

X  =  sqrt((m-l)^2  +  (n-l)^2); 
if  X  <=  r; 

temp(m,n)=besseln(0,a*x) ; 
end; 
end; 

end; 


% 

y (N0:N,N0;N) 
Y(2;N0,N0;N) 
y(2:N0,2:N0) 
Y(N0;N,2:N0) 


=  temp; 

=  flipud(temp) ; 

=  rot90(temp,2) ; 
=  fliplr(temp) ; 


%%%%%For  testing  purposes;  mesh(y) 
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